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Abstract 

Green’s  Function  Analysis  of  Pavement  Deflections 
due  to  Moving  Wheel  Loads 

by  Richard  G.  Thuma 

Chairperson  of  the  Supervisory  Committee:  Professor  Steven  L.  Kramer 

Department  of  Civil  Engineering 


Nondestructive  pavement  deflection  testing  is  commonly  used  to  evaluate  the 
structural  capacity  of  flexible  pavements.  Different  types  of  tests,  characterized  by  their 
loading  mechanism  —  static  or  dynamic,  stationary  or  moving  —  are  used  for  this  purpose. 
For  all  tests,  a  pavement  deflection  model  is  used  to  interpret  the  measured  deflections. 
This  thesis  presents  the  conceptual  development,  validation,  and  parametric  study  of  a 
flexible  pavement  model  for  predicting  deflections  caused  by  a  rolling  wheel. 

The  model  uses  a  Green’s  function  approach  to  compute  deflections.  The  Green’s 
function  computes  the  deflections  at  a  known  distance  from  a  disk  that  is  loaded 
harmonically  at  a  known  frequency.  The  deflections  caused  by  an  impulse  load  are 
computed  using  Fourier  superposition  analysis.  The  rolling  wheel  load  is  modeled  as  a 
sequence  of  impulse  loads  and  the  model  combines  the  responses  from  each  impulse  to 
predict  a  deflection  basin  that  would  be  seen  by  an  observer  moving  with  the  wheel.  This 
process  assumes  that  all  materials  in  the  pavement  structure  behave  linear-elastically. 

The  model  was  validated  by  predicting  deflections  caused  by  falling  weight 
deflectometer  (FWD)  tests.  Predicted  deflections  were  reasonably  close  to  measured 
deflections;  however,  the  pavement  layer  moduli  used  in  the  validation  were 
backcalculated  using  static  rather  than  dynamic  models. 

The  sensitivity  of  the  computed  deflection  basin  to  three  categories  of  input 
parameters  --  algorithm,  material,  and  loading  —  was  investigated  by  a  detailed 
parametric  study.  Algorithm  variables  were  studied  to  determine  the  values  that 


produced  the  most  accurate  solution  in  the  minimum  time.  The  sensitivity  of  the  basin 
shape  to  layer  material  variables  provided  insight  into  which  variables  can  be  potentially 
backcalculated  with  this  model.  Loading  variables  were  studied  to  determine  the  effect 
that  the  wheel  speed,  load,  and  tire  pressure  had  on  the  deflection  basin  shape.  The 
parametric  study  revealed  that  backcalculation  of  layer  moduli  or  thicknesses  may  be 
possible  using  this  model.  However,  backcalculating  both  layer  modulus  and  thickness 
simultaneously  will  require  a  very  precise  definition  of  the  deepest  part  of  the  deflection 
basin.  The  model  may  not  be  practical  for  use  in  a  backcalculation  program  until  the 
computation  time  can  be  reduced. 
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Chapter  1 
INTRODUCTION 


1.1  Background 

Researchers  have  recognized  for  many  years  that  the  structural  capacity  and  life 
expectancy  of  roadways  was  related  to  the  load-deflection  properties  of  the  pavement 
system  materials.  To  determine  these  properties,  two  basic  testing  methods  can  be  used  — 
destructive  or  non-destructive. 

Destructive  testing  requires  removing  material  from  the  roadway  and  then 
performing  laboratory  tests  on  the  samples.  An  obvious  disadvantage  of  this  method  of 
testing  is  the  time  and  cost  associated  with  collecting  samples  and  repairing  the  damaged 
pavement.  However,  there  are  two  very  important  additional  disadvantages.  First, 
samples  cut  from  the  roadway  and  underlying  soils  are  considerably  disturbed  —  this  casts 
doubt  on  the  transferability  of  the  results  of  lab  tests  to  the  field  situation.  Second,  the 
roadway  is  composed  of  several  layers  of  materials  (illustrated  in  Figure  1-1),  each 
composed  of  inherently  heterogeneous  materials.  Destructive  testing  considers  only  a 
small  sample  of  the  materials,  and  does  not  test  the  interaction  between  material  layers. 

Non-Destructive  testing  (NDT),  as  the  name  suggests,  does  not  require  sample 
collection.  Instead  it  indirectly  measures  material  properties  in  situ.  Besides  the  shorter 
test  times  and  relatively  low  cost,  NDT  is  superior  to  destructive  testing  because  it 
measures  a  larger  volume  of  material,  and  tests  the  way  the  individual  layers  interact.  An 
additional  advantage  to  non-destructive  testing  is  that  the  pavement  can  be  evaluated 
using  the  same  load  that  it  is  designed  to  resist  -  i.e.  the  pavement  can  be  loaded  with  a 
truck. 

Several  types  of  non-destructive  tests  have  been  developed  to  measure  the 
deflections  caused  by  an  applied  load.  Static  tests  apply  a  constant  load  to  a  stationary 
area  and  dynamic  tests  apply  a  vibrating  or  impact  load  to  a  stationary  contact  plate. 

Each  test  measures  the  associated  static  or  dynamic  deflections  using  a  variety  of  devices. 
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Travel  Travel 


Figure  1-1.  Typical  three-layered  flexible  pavement  cross  section. 

Recently  a  new  type  of  testing  device  has  been  developed  that  measures  the 
deflections  caused  by  a  rolling  wheel.  This  type  of  device  has  tremendous  potential  for 
determining  the  true  dynamic  properties  of  pavement  structures  because  it  applies  exactly 
the  same  load  that  is  applied  by  highway  vehicles.  Using  this  device,  a  better 
understanding  of  the  deflection  behavior  and  resulting  damage  caused  by  traffic  can  be 
developed. 

All  non-destructive  tests  require  models  to  interpret  the  measured  deflection  data 
because  the  material  properties  cannot  be  measured  directly  (i.e.  deflection  is  not  a 
material  property).  Using  models,  the  layer  material  properties  are  adjusted  until  the 
predicted  deflections  match  the  measured  deflections.  Using  a  process  called 
backcalculation,  the  layer  material  properties  are  adjusted  until  the  deflections  predicted 
by  the  model  match  the  measured  deflections.  The  properties  used  to  calculate  the 
matching  prediction  then  correspond  to  the  properties  of  the  pavement  structure. 

Models  range  in  complexity  from  simple  static  load  formulas  to  intricate  dynamic 
finite  element  computer  programs.  Many  static  and  several  dynamic  models  have  been 
developed  to  interpret  stationary  load  test  deflections.  However,  only  a  few  models  exist 
that  predict  deflections  caused  by  a  moving  load. 
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1.2  Objectives 

The  main  objective  of  this  research  was  to  develop  a  model  to  predict  the 
deflections  resulting  from  a  rolling  wheel  load.  To  accurately  model  pavement  structures, 
the  model  must  be  capable  of  considering  a  multi-layered  system.  Towards  this  end,  a 
Green’s  fimction  approach  was  selected  as  the  method  for  dynamic  analysis.  Also,  the 
model  must  be  capable  of  predicting  the  shape  of  a  constantly-moving  deflection  basin 
because  the  sensors  that  measure  the  deflections  are  attached  to  the  moving  vehicle. 

Because  other  rolling-wheel  models  and  deflection  measurements  were  not 
readily  available,  the  model  was  compared  to  a  different  dynamic  test  —  the  falling 
weight  deflectometer.  This  test  applies  an  impact  load  to  a  stationary  contact  plate.  This 
validation  required  modification  of  the  primary  model  to  allow  for  the  different  loading 
conditions,  but  the  principal  components  remained  identical. 

The  sensitivity  of  the  model  to  the  most  significant  input  parameters  was 
determined  using  a  detailed  parametric  study.  The  changes  in  deflection  basin  shape  due 
to  material  properties  and  loading  conditions  was  investigated.  Also,  the  conceptual 
parameters  that  influence  the  computation  speed  and  accuracy  were  studied  to  determine 
the  values  that  produced  the  best  accuracy  in  the  least  time. 

With  a  completed  model,  backcalculation  can  be  used  to  determine  material 
properties  from  measured  deflections.  However,  developing  a  backcalculation  scheme 
was  outside  the  scope  of  this  research.  Only  the  forward  model  was  created  and 
analyzed. 

1.3  Organization 

The  first  portion  of  this  thesis  describes  the  development  of  non-destructive 
testing  and  the  models  used  to  interpret  the  tests.  The  latter  portion  explains  the 
development  and  validation  of  the  rolling  wheel  model.  The  following  paragraphs 
provide  a  brief  description  of  each  chapter. 
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Chapter  2  describes  the  non-destructive  tests  that  have  historically  been  used  to 
evaluate  pavements.  The  testing  devices  are  grouped  according  to  the  method  by  which 
the  load  is  applied  to  the  pavement:  static,  steady-state  dynamic,  transient  dynamic, 
transient  surface  wave,  and  rolling  wheel.  The  pavement  response  produced  by  each 
loading  scheme  and  the  applications  for  each  test  are  also  described.  The  analytical 
models  used  to  interpret  the  data  measured  by  NDT  devices  are  explained  in  Chapter  3. 
Again,  the  models  are  grouped  by  the  assumed  load  type:  static,  dynamic  stationary  load, 
and  dynamic  rolling  load.  For  each  load  type,  the  simplest  models  are  described  first, 
followed  by  the  more  complicated  models. 

Chapter  4  explains  the  Green’s  fimction  concept  and  describes  how  it  was  used  to 
develop  a  solution  for  the  deflections  of  a  layered  pavement  system  subjected  to 
harmonic  loading.  Chapter  5  provides  a  detailed  conceptual  development  of  an  algorithm 
for  calculating  the  deflection  basin  caused  by  a  rolling  wheel.  The  simplest  case,  a 
stationary  transient  load,  is  considered  first,  followed  by  additional  cases  that  add 
incremental  levels  of  complexity.  Initially,  a  simulated  deflection  response  function  is 
employed  to  easily  illustrate  the  applied  concepts.  Eventually,  the  chapter  describes  how 
the  true  Green’s  function  solution  was  used  with  this  algorithm  in  the  computer  program 
BASIN.  The  program  BASIN  was  verified  by  comparing  the  calculated  deflections  to 
falling- weight  deflectometer  data,  as  described  in  Chapter  6.  The  results  of  a  parametric 
study  of  BASIN  are  presented  in  Chapter  7.  The  study  shows  the  influence  of  the  most 
significant  input  variables  on  the  deflection  basin  shape.  First,  a  “standard”  pavement 
section  was  established,  so  that  all  subsequent  variations  would  have  a  common  basis  for 
comparison.  Next,  each  input  parameter  was  systematically  adjusted  to  illustrate  its 
influence  on  the  calculated  basin.  Finally,  Chapter  8  summarizes  the  findings  of  this 
research  project,  lists  the  conclusions  drawn  from  this  research,  and  makes 
recommendations  for  further  study. 

The  model  developed  in  this  research  was  coded  into  a  computer  program  called 
BASIN.  The  FORTRAN  code  for  BASIN  is  listed  in  Appendix  A. 


Chapter  2 

NONDESTRUCTIVE  PAVEMENT  DEFLECTION  ANALYSIS 

The  first  nondestructive  pavement  testing  device,  the  Travel  Gage,  was  installed 
by  General  Electric  in  1938  on  a  California  highway.  It  consisted  of  holes  drilled  in  the 
pavement  to  different  depths  with  calibrated  rods  inserted  in  the  holes.  As  a  vehicle 
drove  over  the  instrumented  pavement  section,  technicians  observed  and  manually 
recorded  the  deflections  (Moore  et  al.,  1978).  Although  this  device  performed  reasonably 
well,  engineers  soon  realized  that  a  mobile,  non-intrusive  pavement  testing  device  would 
be  much  more  useful  for  assessing  pavement  and  subgrade  stiffiiess. 

As  numerous  technological  and  theoretical  advances  occurred  over  the  next  50 
years,  engineers  developed  many  different  devices  for  determining  the  stiffness  of  a 
pavement  system.  These  nondestructive  devices  can  be  grouped  into  five  broad 
categories,  distinguished  by  their  different  loading  schemes. 

1)  Static 

2)  Steady-State  Dynamic 

3)  Transient  Dynamic 

4)  Transient  Surface  Wave  Analysis 

5)  Moving  Wheel 

This  chapter  describes  the  general  concepts,  loading  responses,  testing  equipment,  and 
applications  for  each  category  of  pavement  testing  device. 

2.1  Static  Load 

2.1.1  General 

Pavement  evaluation  systems  using  stationary  or  very  slowly  moving  loads  fall 
into  the  category  of  static  deflection  devices.  Measurements  are  taken  before  and  after 
loading,  or  at  points  radially  outward  from  a  loaded  plate.  The  primary  problem  with 
these  systems  is  establishing  a  reference  point  from  which  to  measure  the  deflection  --  the 
point  must  be  close  enough  to  be  convenient,  yet  far  enough  away  that  it  is  not  affected 
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by  the  load.  Because  of  this  problem,  the  deflection  data  obtained  by  static  devices  is 
often  questionable.  However,  because  so  much  data  has  been  collected  with  these 
devices,  the  developed  correlations  are  considered  fairly  reliable  (Moore  et  al.,  1978) 

The  first  static  deflection  measurement  device,  the  GE  Travel  Gage,  was  intrusive, 
using  rods  inserted  into  drilled  holes  of  different  depths  to  measure  the  deflection  of 
pavement  and  subgrade  layers.  To  improve  data  acquisition  speed,  the  rods  were 
outfitted  with  linear  variable  differential  transformers  (LVDTs  -  an  electronic  dial  gauge). 
This  device  was  extensively  used  to  measure  pavement  deflections  during  the  Western 
Association  of  State  Highway  Officials  (WASHO)  Road  Test  conducted  in  southeastern 
Idaho  in  1952.  However,  the  Travel  Gage  still  proved  too  slow  and  expensive  for 
widespread  use  (Moore  et  al.,  1978). 

Based  on  observations  at  the  WASHO  Road  Test,  planners  realized  a  mobile 
pavement  testing  device  would  be  needed  to  evaluate  pavement  at  specific  locations 
while  the  test  progressed.  Fortunately,  A.C.  Benkelman  was  developing  a  simple  lever- 
arm  device  to  measure  static  deflections  of  flexible  pavements.  He  released  this  device, 
dubbed  the  Benkelman  Beam,  in  1952.  As  early  as  1954,  the  California  Division  of 
Highways  began  widespread  use  of  the  Benkelman  Beam  (Moore  et  al.,  1978).  Many 
other  agencies  followed  California’s  lead,  while  others  chose  to  improve  upon  the 
manually  operated,  tedious  instrument.  These  versions,  using  the  same  concept,  included 
the  Traveling  Deflectometer  and  the  LaCroix  -  LiC.P.C.  (Laboratoire  Central  des  Fonts  et 
Chaussees)  Deflectograph,  which  greatly  improved  the  data  collection  speed  (Moore  et 
al.,  1978). 

2.1.2  Pavement  Deflection  Response 

Static  loading  devices  create  a  static  deflection  basin,  which  has  a  magnitude  and 
shape  that  is  dependent  on  the  stiffness  of  the  pavement  system.  Deflection  basins 
produced  by  static  testing  devices  are  relatively  large  (up  to  10.2  mm  (400  mils)) 
compared  to  other  testing  devices  because  the  loads  are  large  and  applied  for  a  long  time 
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(several  minutes).  (Hoffman  and  Thomson,  1981).  Very  slow  moving  loads  (less  than 
3.2  km/hr  (2  mph))  create  a  moving  deflection  basin  that,  for  all  practical  purposes,  is 
identical  to  a  stationary  basin.  Thus,  static  devices  have  two  potential  methods  for 
determining  the  shape  of  the  deflection  basin. 

1)  A  stationary  load  can  be  applied  and  deflections  measured  directly  under 
and  at  points  radially  outward  from  the  loading  plate. 

2)  Deflections  can  be  measured  periodically  at  one  stationary  point,  while  the 
load  is  slowly  rolled  away. 

The  deflections  obtained  from  either  of  these  methods  are  plotted  as  function  of 
distance  from  the  load,  as  illustrated  in  Figure  2-1 .  Rotating  this  deflection  shape  around 
the  center  of  the  load  creates  the  entire  static  deflection  basin.  Once  the  shape  of  the 
basin  is  defined,  parameter  correlations  or  backcalculation  techniques  are  used  to 
determine  the  pavement  system  properties. 


Figure  2-1.  Deflection  basin  and  individual  deflections  produced  by  a  loaded  plate  as  measured  by  seven 
deflection  sensors. 


Correlations  directly  relate  one  measurement  or  index  to  a  pavement  property.  In 
pavement  analysis,  two  parameters  are  commonly  used  for  correlations: 
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D(r)  -  the  magnitude  of  the  deflection,  D,  at  a  distance,  r,  from  the 
center  of  the  loading  plate.  Many  correlations  between  deflections  and  moduli 
have  been  developed,  and  they  generally  follow  the  trend  that  deflections 
under  and  nearby  the  load  plate  correlate  best  to  AC  layer  modulus, 
deflections  in  the  middle  of  the  basin  (between  300  and  900  mm  (1 1 .8  and 
35.4  in)  from  the  center  of  the  load)  correlate  best  to  the  base  course  modulus, 
and  deflections  at  the  edge  of  the  basin  relate  best  to  the  subgrade  modulus. 

Area  Parameter  —  a  normalized  index  equal  to  the  sum  of  the  areas  of  the 
trapezoids  deflned  by  the  deflection  measurements  divided  by  the  deflection 
directly  under  the  center  of  the  loading  plate.  For  example,  the  area  parameter 
from  Figure  2-1  would  be  given  by: 

A  =  —(Do  +  2D,  +  2D2  +  2D3  +  2D4  +  2D5  +  Dg) 

^0 

Several  correlations  have  been  made  between  the  area  parameter  and 
pavement  layer  stiffiiesses. 

Backcalculation  techniques  are  iterative,  require  multiple  calculations,  and  use 
computer  programs  to  quickly  solve  complicated  mathematics.  First,  a  predicted 
deflection  basin  is  calculated  using  a  model  and  compared  to  the  measured  deflection 
basin.  Next,  certain  material  parameters  are  adjusted,  and  Anally  the  predicted  deflection 
basin  is  recalculated,  beginning  the  next  iteration.  Iterations  continue  until  the  predicted 
deflection  basin  matches  the  measured  deflection  basin,  within  a  tolerable  error.  The 
values  of  material  parameters  that  were  used  to  produce  the  matching  deflection  basin 
correspond  to  the  actual  payement  parameters. 

2.1.3  Testing  Devices 

ASTM  D4695-87  (Standard  Guide  for  General  Pavement  Deflection 
Measurements)  provides  guidance  for  using  static  deflection  devices. 
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2. 1 .3 . 1  Benkelman  Beam 

A  long  beam  carefully  placed  between  the  dual  tires  of  a  loaded  truck  (usually 
with  80  kN  (1 8,000  lb)  rear  axle  weight)  measures  the  deflected  pavement  position.  As 
the  truck  slowly  drives  away  and  the  pavement  rebounds,  the  beam  rotates  around  a 
stationary  pivot  point.  The  new  position  represents  the  undeflected  pavement  elevation 
and  is  compared  to  the  deflected  elevation.  An  observer  reads  the  measurements  from  a 
dial  gauge  that  is  accurate  to  25.4  pm  (0.001  in)  and  then  calculates  the  deflection. 

Figure  2-2  illustrates  the  basic  components  of  the  Benkelman  Beam. 

The  Benkelman  Beam  is  simple  and  inexpensive,  but  relatively  slow.  The 
instrument  alone  costs  only  about  $1,000  and  an  experienced  crew  can  make  between  300 
and  400  measurements  in  a  day.  Of  course,  traffic  must  be  carefiilly  controlled  around 
the  testing  crew  to  ensure  their  safety. 


Figure  2-2.  Schematic  representation  of  the  Benkelman  Beam  (after  Meier,  1995). 


2. 1 .3 .2  Traveling  Deflectometer 

The  Traveling  Deflectometer  is  essentially  an  automated  version  of  the 
Benkelman  Beam  and  was  used  by  the  California  Division  of  Highways.  A  tractor-trailer 
outfitted  with  two  lever  arms  simultaneously  measures  the  deflection  between  each  set  of 
dual  tires  on  its  80  kN  (18,000  lb)  rear  axle.  Undeflected  measurements  are  recorded  by 
separate  probes  located  beyond  the  loaded  axle’s  influence  zone.  Measurements  are 
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recorded  every  6.1  m  (20  ft)  as  the  vehicle  rolls  forward  at  a  creep  speed  of  0.8  km /hr 
(0.5  mph). 

With  this  device,  1,500  to  2,000  measurements  are  possible  in  a  single  day 
(Moore  et  al.,  1978).  Also,  safety  was  greatly  improved  with  this  device  because  the 
operators  ride  on  the  vehicle  during  the  test. 

2. 1 .3 .3  LaCroix  -  L  C.P.  C.  Deflectograph 

The  LaCroix  -  L.C.P.C.  Deflectograph  is  very  similar  to  the  Traveling 
Deflectometer  —  it  also  is  essentially  an  automated  Benkelman  Beam.  It  uses  the  1 16  kN 
(26,000  lb)  re^  axle  of  a  truck  to  load  the  pavement,  and  measures  deflections  with  a 
unique  moving  beam  system.  As  the  truck  rolls  slowly  forward  at  1 .8  km/hr  (1.1  mph), 
twin  measurement  beams,  each  analogous  to  a  Benkelman  Beam,  are  lowered  to  rest  on 
the  pavement  surface.  The  frame  supporting  the  beams  slides  slowly  backwards  to  keep 
the  probe  tips  in  exactly  the  same  place  while  the  wheels  roll  forward.  The  probe  tips 
start  about  one  meter  in  .front  of  the  rear  wheels  and  stay  on  the  pavement  imtil  they  reach 
a  maximum  deflection,  slightly  behind  the  rear  wheels.  Thus,  entire  deflection  basins  can 
be  measured  and  automatically  recorded  by  the  Deflectograph.  When  a  test  is  complete, 
the  beams  quickly  slide  forward  to  start  the  testing  cycle  again  (Prandi,  1967). 

The  LaCroix  -  L.C.P.C  Deflectograph  is  significantly  faster  than  the  Benkelman 
Beam  or  the  Traveling  Deflectometer.  Prandi  (1967)  reported  that  7,000  to  8,000 
measurements  were  possible  in  a  single  day.  However,  he  added  a  caveat  to  this  report 
that  speaks  to  the  accuracy  of  the  measurements,  “...pavement  deformability  is  better 
expressed  through  analysis  of  numerous  deflection  measurements,  even  not  very  accurate, 
than  by  that  of  a  few  measurements,  even  highly  accurate.” 

2. 1.3. 4  Plate  Bearing  Test 

In  the  plate  bearing  test,  a  hydraulic  jack  reacts  against  a  stationary  frame  (usually 
a  heavy  tractor-trailer)  to  load  a  steel  plate.  This  method  is  very  slow  (2-3  hours  per  test) 
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and  used  primarily  for  asphalt  overlay  thickness  design  for  existing  pavement 
structures.  The  US  Navy  procedure  for  this  test  uses  two  separate  load  plates,  203  and 
760  mm  (8  and  30  in)  diameter,  and  loads  are  applied  until  a  5.1  mm  (0.2  in)  deflection  is 
observed  directly  under  the  load  plate  (Moore  et  al.,  1978). 

2.1.4  Application 

Static  deflection  devices  showed  promise  and  filled  a  valid  role  during  their  time, 
but  now  have  been  outdated  by  technological  advances.  However,  the  data  relationships 
developed  for  static  devices  still  are  used  today  for  correlating  steady-state  and  impulsive 
loading  testing  to  pavement  overly  design. 

Even  though  static  devices  were  widely  used  by  highway  officials,  engineers 
realized  that  static  loading  devices  had  several  problematic  features.  First,  deteriorated, 
loosely  bound  pavement  surfaces  are  difficult,  if  not  impossible,  to  test  with  static  devices 
because  the  probe  tip  can  move  erratically  if  it  rests  on  a  loose  stone.  Second,  a 
convenient,  undeflected  reference  point  for  deflection  measurements  is  difficult  to  find. 
Stiffer  pavements  with  bound  base  course  materials  have  wider  basins  and  smaller 
amplitude  deflections,  thus,  the  Deflectograph  and  Traveling  Deflectometer  may  not 
record  reference  points  that  are  truly  undeflected.  Third,  the  static  loading  scheme  does 
not  incorporate  the  dynamic  effects  of  a  rolling  wheel  load. 

2.2  Steady  State  Dynamic  Load 

The  next  generation  of  pavement  testing  devices,  steady-state  dynamic  load, 
overcame  some  of  the  limitations  of  static  testing  devices. 

2.2.1  General 

All  steady-state  deflection  devices  cyclically  load  and  unload  the  pavement  to 
induce  a  steady-state  sinusoidal  vibration  in  the  pavement.  The  dyneimic  force  generator 
applies  a  load  time  history  similar  to  that  shown  in  Figure  2-3  by  superimposing  the 


12 


dynamic  load  on  the  static  load  produced  by  the  weight  of  the  equipment.  Note  that 
the  dynamic  force  (peak-to-peak)  must  be  less  than  twice  the  static  weight  of  the 
equipment  to  ensure  continuous  contact  with  the  pavement  surface. 

Steady-state  devices  represent  a  significant  improvement  over  static  instruments 
because  they  overcome  the  problem  of  finding  a  stationary  reference  point  for  measuring 
deflections.  Rather  than  measuring  the  absolute  deflection  due  to  the  dynamic  load,  the 
device  measures  the  amplitude  of  the  relative  deflection  produced  by  the  dynamic  load. 

In  other  words,  an  inertial  reference  system  is  used  and  the  amplitude  of  the  deflection 
(peak-to-peak  deflection  value)  can  be  compared  directly  to  the  amplitude  of  the  dynamic 
force  (peak-to-peak  force  value).  Therefore,  for  a  given  dynamic  force  amplitude,  a 
smaller  deflection  amplitude  indicates  a  stiffer  pavement. 


Time  - > 


Figure  2-3.  Typical  output  of  force  generator  used  for  steady-state  dynamic  testing. 

Improvements  in  deflection  measurement  systems  also  make  steady-state  devices 
a  practical  alternative  to  static  systems.  With  steady-state  devices,  deflections  are 
measured  with  electronic  displacement,  velocity,  or  acceleration  transducers,  rather  than 
dial  gauges  or  LVDTs.  These  inertial  motion  sensors  are  essentially  single  degree  of 
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freedom  structures.  An  internal  mass  moves  within  a  casing  along  a  particular  axis 
(vertical,  for  pavement  deflection  devices).  The  internal  mass  displacement,  velocity,  or 
acceleration  relative  to  the  external  case  generates  an  output  signal  that  is  processed  and 
recorded  by  a  computer  system  (Graves  and  Dmevich,  1993). 

Most  steady-state  devices  use  velocity  transducers  (commonly  called  geophones) 
and  an  automatic,  electronic  integrator  to  convert  the  velocity  time  history  to  a 
displacement  time  history.  Velocity  transducers  with  integrators,  rather  than 
displacement  transducers,  are  commonly  used  because  the  output  characteristics  of  most 
displacement  transducers  are  heavily  frequency  dependent,  whereas  velocity  outputs  are 
frequency  independent  above  a  minimum  frequency  (Moore  et  al.,  1978). 

2.2.2  Pavement  Deflection  Response 

In  static  loading  systems,  when  the  load  is  applied,  the  pavement  deflects 
proportionally  to  the  applied  load;  when  the  load  is  removed,  the  pavement  returns  to  the 
same  position  (assuming  the  load  is  within  the  elastic  range).  The  same  principles  apply 
to  dynamic  loading.  The  oscillating  load  creates  a  maximum  deflection  when  the 
maximum  load  is  applied,  and  the  pavement  rebounds  as  the  load  amplitude  reaches  a 
minimum.  Therefore,  the  amplitude  of  the  dynamic  deflection  is  proportional  to  the 
dynamic  force.  In  effect,  the  load  on  the  pavement  creates  a  vibration  in  the  pavement 
oscillating  at  the  same  frequency  that  the  load  is  applied.  The  amplitude  of  the  deflection 
decreases  as  the  distance  between  the  load  and  the  measurement  location  increases 
because  of  material  and  radiation  damping. 

The  deflection  basins  measured  in  a  steady-state  dynamic  test  are  similar  to  the 
deflections  created  by  a  static  load.  However,  dynamic  deflections  are  different  because 
the  dynamic  stiffness  of  the  pavement  system  changes  with  the  driving  frequency.  The 
dynamic  stiffness  has  been  found  to  increase  as  the  driving  frequency  increases.  Thus, 
smaller  deflections  are  observed  during  higher  frequency  tests  (Moore  et  al.,  1978). 
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2.2.3  Tggting  E^quipment 

ASTM  4695-87  (Standard  Guide  for  General  Pavement  Deflection 
Measurements)  describes  the  general  procedures  for  testing  pavements  with  steady-state 
dynamic  devices. 

2.2.3. 1  Dynaflect 

The  Dynaflect  consists  of  several  deflection  sensors  and  a  motor  that  drives  two 
coimter-rotating  masses,  mounted  on  a  two- wheel  trailer.  The  resulting  dynamic  load  is 
transferred  to  the  pavement  through  two  101  mm  (4  in)  wide,  406  mm  (16  in)  outer- 
diameter  rubber-coated  steel  wheels  spaced  504  mm  (20  in)  center  to  center.  The 
resulting  pavement  contact  area  is  less  than  2581  mm  (4  in  ).  Figure  2-4  illustrates  a 
general  schematic  of  the  device  and  the  typical  loading  configuration. 

When  the  device  stops  at  the  test  point,  the  sensors  are  lowered  to  the  pavement 
on  a  beam  aligned  with  the  axis  of  symmetry  which  passes  between  the  wheels. 
Geophones  spaced  304.8  mm  (12  in)  apart  record  the  deflections  induced  by  the  dynamic 
force  generator.  The  machine’s  static  load  of  8  kN  (1,800  lb)  is  offset  by  the  peak-to- 
peak  dynamic  force  of  4.4  kN  (1,000  lb)  applied  at  a  fixed  frequency  of  8  Hz. 

The  driver  of  the  towing  vehicle  completely  controls  the  Dynaflect  via  remote 
console.  The  entire  operation  at  each  test  site  obtains  5  measurements  within  about  two 
minutes  (Moore  et  al.,  1978).  Because  the  device  is  stationary  during  testing,  traffic  must 
be  controlled  in  the  vicinity  of  the  machine. 

2.2.3 .2  Road  Rater 

Several  models  of  this  instrument  are  available,  but  they  all  share  certain  common 
features.  First,  a  vibratory  load  is  produced  by  a  steel  mass  and  hydraulically  actuated 
vibrator.  The  pavement  is  loaded  through  two  102  x  1788  mm  (4x7  in)  steel  plates 
spaced  267  mm  (10.5  in)  center  to  center,  for  a  total  contact  area  of  36,129  mm  (56  in  ). 
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b)  Plan  View 
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Figure  2-4.  a)  Schematic  Illustration  of  the  Dynaflect.  b)  Plan  view  of  the  loading  and  geophone 
configuration. 


Normally  four  geophones,  spaced  at  305  mm  (12  in)  along  the  S5mimetry  axis 
passing  between  the  load  plates,  measure  the  deflections.  Like  the  Dynaflect,  these  are 
lowered  to  the  pavement  at  each  test  location.  Depending  on  the  model,  the  Road  Rater 
can  apply  loads  ranging  between  8  and  22  kN  (1,800  and  5,000  lb)  at  five  different 
frequencies,  ranging  between  10  and  40  Hz  (Moore  et  al.,  1978). 

Tests  are  very  fast,  requiring  less  than  a  minute  for  one  deflection  measurement  at 
a  preset  driving  frequency.  Like  the  Dynaflect,  the  Road  Rater  is  operated  by  the  driver 
of  the  towing  vehicle. 


2.2.4  Application 

Steady-state  deflection  testing  devices  overcome  some  of  the  problems  with  static 
devices.  The  most  important  advance  is  the  use  of  an  inertial  reference  system  for 
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measuring  deflections.  As  described  in  Section  2.2.1,  very  sensitive  and  accurate 
accelerometers  or  velocity  transducers  indirectly  measure  deflections  of  the  pavement 
surface.  Because  a  stationary  local  reference  is  not  needed,  the  devices  are  much  smaller, 
lighter,  and  consequently  more  mobile.  All  this  results  in  cost  savings  for  highway 
departments,  allowing  a  smaller  crew  (1  person)  to  accomplish  more  measurements  in 
less  time. 

Even  though  steady-state  testing  devices  are  a  significant  improvement  over  static 
load  devices,  there  are  still  problems  with  the  test.  First,  moving  traffic  does  not  load  the 
pavement  sinusoidally  like  steady-state  devices.  Second,  the  equipment  must  be 
stationary  during  the  test,  which  requires  routing  traffic  around  the  equipment.  Third, 
tests  are  only  accomplished  at  finite  locations,  usually  about  every  few  himdred  meters. 
Therefore,  a  weak  section  could  go  undetected  until  significant  surface  distress  was 
recognized. 

2.3  Transient  (Impulse)  Dynamic  Load 

The  next  generation  of  pavement  testing  devices,  falling  weight  deflectometers, 
attempt  to  overcome  the  single-frequency  loading  problem  associated  with  steady-state 
devices  by  loading  the  pavement  with  a  single  impulse  force. 

2.3.1  General 

All  impact  load  devices,  commonly  referred  to  as  falling  weight  deflectometers 
(FWDs),  produce  a  transient  impulse  load  by  dropping  a  known  weight  from  a  known 
height  onto  a  contact  plate  of  known  area  and  stiffiiess.  Typically,  the  peak  load  and  the 
induced  vertical  deflections  are  measured  at  points  directly  under  the  load  and  at  different 
distances  from  the  load. 

FWD  loading  systems  are  comprised  of  a  force-generating  device  (falling  weight), 
guide  system,  loading  plates,  and  load  cell.  The  falling  weight  is  mechanically  raised  to 
one  of  several  controlled  heights  and  released.  The  guides  ensure  that  the  weight  falls 
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perpendicular  to  the  pavement  surface  with  negligible  friction,  and  that  the  weight 
strikes  the  plate  in  a  controlled,  repeatable  manner.  Load  plates  help  to  xmiformly 
distribute  the  load  over  a  305  to  457  mm  (12  to  1 8  in)  diameter  area.  They  are  generally 
made  of  steel  with  a  rubber  coating  on  the  bottom,  which  helps  to  condition  and  extend 
the  transient  load  shape.  Additionally,  the  plates  have  a  small  opening  in  their  centers  to 
allow  for  a  deflection  measuring  device.  The  force  pulse  duration  lasts  between  20  and 
60  milliseconds  (ms)  and  can  be  approximately  described  by  a  haversine  or  half-sine 
wave.  This  load  shape  more  closely  resembles  the  action  of  real  traffic  than  either  the 
static  or  steady-state  devices  (Barksdale,  1971). 

Deflection  measurement  sensors  are  generally  mounted  on  a  beam  extending  from 
the  load  plate  radially  outward  along  the  travel  direction.  Each  particular  testing  model 
may  have  different  numbers  and  spacing  of  sensors.  With  some  models,  the  sensor 
locations  may  be  adjusted,  depending  on  the  testing  application.  The  sensors  are 
cormected  to  an  onboard  computer  system  that  records  the  peak  vertical  deflection  at  each 
sensor  location.  Several  types  of  sensors  are  used,  including  seismometers,  velocity 
transducers,  and  accelerometers. 

2.3.2  Pavement  Deflection  Response 

When  the  falling  weight  strikes  the  load  plate,  the  surface  of  the  pavement 
deflects  by  an  amplitude  dependent  on  the  stiffness  of  the  entire  pavement  system.  While 
the  deflection  is  occurring  directly  under  the  load  plate,  a  deflection  wave  front  is 
emanating  radially  away  from  the  load  plate.  As  a  result  of  radiation  and  material 
damping,  the  amplitude  of  the  pavement  deflection  decreases  as  the  wave  moves  farther 
from  the  load  plate.  For  a  given  impulse  force,  a  stiffer  pavement  will  deflect  less  than  a 
softer  pavement.  Also,  the  speed  at  which  the  wave  propagates  away  from  the  load  plate 
is  a  function  of  the  material  stiffness  —  waves  travel  faster  in  stiffer  systems. 

As  described  in  the  previous  section,  FWDs  have  sensors  that  record  the 
deflection  at  several  points  on  a  line  extending  outward  from  the  load  plate.  These 
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sensors  capture  a  time  history  of  deflection  at  each  point;  however,  typically  only  the 
peak  deflection  is  used  for  analysis.  Plotted  together,  these  peaks  form  a  “static” 
deflection  basin.  Close  analysis  of  the  time  history  of  deflections  shows  that  the  peak 
deflections  do  not  occur  simultaneously  at  all  points.  Rather,  the  peak  occurs  later  at 
each  successive  sensor.  Thus,  the  “static”  basin  never  truly  exists  during  an  FWD  test. 
This  concept  is  more  fully  illustrated  in  Chapter  5. 

Once  the  “static”  deflection  basin  is  known,  the  same  procedures  used  for  static 
and  steady-state  devices  are  used  to  analyze  the  data  and  determine  layer  moduli.  That  is, 
correlations  to  basin  measurements  or  backcalculation  techniques  are  used,  as  described 
in  Section  2.1.2. 

2.3.3  Testing  Equipment 

ASTM  4694-87  (Standard  Test  Method  for  Deflections  with  a  Falling  Weight 
Type  Impulse  Load  Device)  describes  the  general  testing  procedures,  apparatus,  and  data 
acquisition  systems  applicable  to  transient  dynamic  load  testing  devices. 

2.3.3. 1  Dynatest  8000  Falling  Weight  Deflectometer. 

Introduced  in  the  early  80’ s,  this  FWD  device  has  the  longest  service  record  in  the 
United  States.  It  is  trailer  mounted  and  remotely  controlled  by  the  driver  of  the  towing 
vehicle.  The  impulse  load  is  delivered  by  dropping  a  hydraulically  raised  weight  onto  a 
buffer/loading  plate  system.  The  loading  plate  is  300  mm  (11.8  in)  in  diameter  and  is 
comprised  of  a  single  piece  of  steel  with  a  rubber  pad  on  the  bottom  to  help  distribute  the 
load  more  evenly.  Figure  2-5  illustrates  the  loading  system. 

Deflections  are  measured  using  7  geophones  located  on  a  raise-lower  bar 
extending  towards  the  towing  vehicle  from  the  load  plate.  An  onboard  computer  records 
the  peak  deflection  or  the  entire  time  history  of  deflection  for  each  sensor.  The  load  pulse 
duration  is  25  to  30  ms,  and  the  load  magnitude  can  range  between  7  and  125  kN  (1,574 
and  28,103  lb),  depending  on  the  drop  height. 
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Up  to  eight  individual  drops  at  4  different  heights  can  be  used  in  a  single  test 
program.  Like  the  steady-state  devices,  the  Dynatest  is  automated  and  requires  less  than 
two  minutes  at  each  test  site. 

2. 3.3.2  KUAB  2M- Falling  Weight  Deflectometer  (Model  8714) 

The  KUAB  2M,  manufactured  in  Sweden,  gives  results  similar  to  the  Dynatest, 
but  loads  and  measures  the  pavement  deflections  in  a  very  different  way.  The  KUAB 
uses  a  two-mass  system,  consisting  of  a  falling  mass  and  a  second  mass-buffer  system  to 
create  the  impulse  load,  as  illustrated  in  Figure  2-5.  The  two-mass  system  creates  a  more 
reproducible  load  shape,  and  also  approximately  doubles  the  pulse  duration  (60  -  70  ms) 
(Crovetti  et  al.,  1989).  The  load  plate  is  a  300  mm  (1 1.8  in)  diameter  steel  plate 
segmented  into  four  quadrants,  each  with  a  rubber  pad.  The  four  segments  are  equally 
loaded,  but  the  segmenting  allows  the  plate  to  better  conform  to  pavement  surface 
irregularities. 


Figure  2-5.  Schematic  illustration  of  the  loading  mechanism  used  by  two  falling  weight  deflectometers 
(after  Crovetti  et  al.,  1995). 
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Rather  than  using  geophones,  the  KUAB  is  equipped  with  seven  seismometers 
(deflection  transducers)  to  measure  the  pavement  displacements.  The  sensors  are 
mounted  on  a  raise-lower  bar  that  extends  radially  away  from  the  towing  vehicle. 

Like  the  Dynatest,  the  KUAB  2M  is  fully  automated,  requiring  just  a  few  minutes 
at  each  test  site.  Unlike  the  Dynatest,  a  virtually  unlimited  number  of  drops  can  be 
programmed  for  each  height  in  a  testing  sequence.  Also,  the  testing  device  is  completely 
enclosed  by  a  locking  metal  housing  for  protection  from  weather  and  road  grime  during 
transportation. 

2.3.4  Application 

Falling  weight  deflectometers  are  widely  used  by  state  highway  administrations 
for  pavement  evaluation  and  overlay  design.  Often,  the  FWD  testing  program  is  used  in 
conjunction  with  a  backcalculation  program  like  EVERCALC,  ILLI-PAVE ,  BISAR,  or 
WESDEF.  These  programs  adjust  layer  stiffnesses  until  the  predicted  deflection  basin 
matches  the  measured  “static”  deflection  basin.  Maintenance  decisions  and  rehabilitation 
design  are  based  on  the  layer  stiffiiesses  determined  by  the  FWD  testing  program. 

Some  researchers  have  developed  methods  for  determining  layer  thickness  along 
with  stiffhess  (Sivaneswaran  et  al.,  1993).  Additionally,  more  work  is  xmderway  that  uses 
the  entire  time  history  of  deflection  for  pavement  analysis  (Ketchum,  1993;  Zaghloul  et 
al.,  1994). 

Although  deflection  analysis  programs  are  now  well-developed  and  produce 
accurate,  repeatable  results,  there  are  still  two  principle  limitations  to  FWD  testing.  First, 
the  FWD  must  be  completely  stationary  to  perform  each  test.  Therefore  traffic  control  is 
needed  to  keep  traffic  flowing  aroimd  the  machine.  Second,  although  the  FWD  load 
pulse  is  very  similar  to  a  rolling  wheel  passing  over  a  stationary  point,  it  is  not  exactly  the 
same.  Coupled  with  analysis  programs  that  only  use  peak  “static”  deflection  basins,  the 
true  dynamic  response  of  the  pavement  is  not  well  understood.  In  spite  of  the 
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disadvantages  of  the  FWD,  it  remains  far  superior  to  earlier  testing  devices.  As 
backcalculation  models  become  more  refined,  results  will  only  improve. 

2.4  Transient  Surface  Wave  Analysis 

2.4.1  General 

A  fourth  type  of  nondestructive  testing  called  Spectral  Analysis  of  Surface  Waves 
(SASW)  can  be  used  to  determine  the  pavement  system  layer  modulus  and  thickness. 
Like  falling  weight  deflectometers,  SASW  uses  the  results  of  deflection  time  histories  at 
discrete  points  due  to  an  impulse  load.  However,  instead  of  using  deflection  magnitudes, 
the  SASW  method  uses  the  phase  dijference  in  surface  wave  records  to  determine  the 
wave  propagation  velocity  at  various  depths.  When  the  velocity,  V^,  density,  p,  and 
Poisson’s  ratio,  v,  are  known,  the  material  modulus  (a.k.a.  stiffness),  E,  can  be  directly 
calculated  through  the  following  basic  relationships, 

E  =  2G(l  +  v)=2pV/(l+y) 
where  G  is  equal  to  the  shear  modulus. 

Figure  2-5  schematically  illustrates  the  basic  components  used  in  SASW  testing. 
The  impulse  source,  ranging  from  a  tack  hammer  to  a  67  N  (15  lb)  sledge  hammer 
(depending  on  the  distance  to  the  geophones),  impacts  the  pavement  and  generates  a 
group  of  Rayleigh  waves  (surface  waves)  emanating  from  the  source  at  various 
frequencies.  The  two  vertically  oriented  receivers  record  the  motion  of  the  pavement  as 
the  wave  energy  moves  past  them.  The  spectral  analyzer  quickly  converts  the  time 
history  output  from  the  receivers  into  a  frequency  domain  response  using  a  fast-Fourier 
transformation  (FFT).  By  analyzing  the  phase  difference  between  the  two  receivers  at 
each  frequency,  the  variation  of  Rayleigh  wave  velocity  with  wavelength  is  determined. 

A  plot  of  the  surface  wave  velocity  versus  the  wavelength  is  called  a  dispersion  curve. 
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The  dispersion  curve  is  the  “raw  data”  carried  out  of  the  field  for  in-house  data 
reduction  (Nazarian  and  Stokoe,  1984).  . 

Dispersion  curves  are  then  converted  to  meaningful  layer  stiffness  and  thickness 
information  by  a  process  knovm  as  inversion.  The  earliest  and  least  sophisticated 
inversion  processes  simply  assumed  that  the  shear  wave  velocity  was  equal  to  1.1  times 
the  Rayleigh  wave  velocity  at  the  depth  corresponding  to  one-half  to  one-third  the 
wavelength  (Nazarian  and  Stokoe,  1984).  In  recent  years,  more  sophisticated  inversion 
processes  have  been  developed  that  model  the  pavement  as  a  laminate  of  many  individual 
sublayers,  each  with  (potentially)  imique  properties.  Through  an  iterative  process,  a  shear 
wave  velocity  (or  stiffness)  for  each  sublayer  is  assumed,  and  a  theoretical  dispersion 
curve  is  calculated.  The  theoretical  curve  and  the  measured  curve  are  compared,  and  then 
the  shear  wave  velocity  profile  is  adjusted  until  the  two  curves  match,  within  a  tolerable 
error. 


Spectral  Analyzer 


Figure  2-6.  General  configuration  of  SASW  test  (after  Nazarian  and  Stokoe,  1989). 
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The  final  shear  wave  velocity  profile  is  converted  to  a  stiffness  profile  by 
assuming  a  Poisson’s  ratio  and  mass  density  for  each  layer,  and  using  the  formula: 

i:=2G(l+v)=2pF/(H-v) 

The  stiffhess  profile  is  next  examined  to  determine  layer  thicknesses.  The 
boundary  between  separate  layers  is  indicated  by  a  sharp  contrast  in  stiffness.  Thus,  if 
two  adjacent  sublayers  have  very  similar  stiffnesses,  they  can  be  considered  as  one  layer 
(Nazarian  and  Stokoe,  1989). 

2.4.2  Pavement  Deflection  Response 

As  established  in  the  previous  section,  SASW  tests  derive  all  of  their  information 
from  the  propagation  of  surface  waves  over  a  broad  range  of  frequencies.  If  the  SASW 
test  were  conducted  on  a  homogeneous,  isotropic,  elastic  halfspace,  only  one  frequency 
would  be  needed  for  testing,  because  the  surface  wave  velocity  would  be  constant  for  all 
frequencies.  However,  because  pavement  layer  properties  vary  with  depth,  the  surface 
wave  is  dispersive;  i.e.  its  velocity  changes  with  frequency.  Thus,  to  determine  the 
variation  of  stiffness  from  the  pavement  surface  down  through  the  subgrade,  tests  must  be 
conducted  using  very  high  (tack  hammer)  to  very  low  (sledge  hammer)  impulse 
frequencies. 

In  theory,  all  frequencies  should  be  recorded  with  a  single  impact.  However,  in 
practice,  certain  receiver  spacings  and  impulse  loadings  work  best  to  develop  a  quality 
profile.  Nazarian  and  Stokoe  (1989)  describe  a  common  receiver  midpoint  (CMRP)  array 
that  is  most  appropriate  for  analyzing  pavement  layers.  Forward  and  reverse  profiles  are 
conducted  with  each  configuration  to  reduce  data  scatter.  Figure  2-7  illustrates  the  source 
and  receiver  arrangement  for  the  CMRP  array. 

2.4.3  Testing  Equipment 

The  field  equipment  required  for  SASW  testing  basically  consists  of  three 
components:  a  source,  two  (or  more)  receivers,  and  a  recording/analyzing  device. 
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Sources  used  in  S  ASW  testing  must  be  capable  of  generating  surface  waves 
containing  a  wide  range  of  frequencies.  Because  small  hammers  generate  predominantly 
high  frequency  waves  and  large  hammers  generate  predominantly  low  frequency  waves,  a 
combination  of  hammers  is  often  used  for  the  entire  series  of  tests  at  one  site.  Extremely 
high  frequency  waves  (over  20  kHz)  must  be  generated  to  examine  the  upper  few 
centimeters  of  pavement;  pulsating  crystals  have  been  successfully  used  for  this  purpose 
(Nazarian  and  Stokoe,  1989).  Very  small  impulse  sources  can  be  used  when  the  sensors 
are  close  together  because  the  deflection  amplitude  is  not  used  for  S ASW  testing,  only 
the  phase  lag  between  receivers  is  important. 

Receivers  are  selected  depending  on  the  frequency  of  interest.  Accelerometers  are 
better  at  capturing  higher  frequency  waves  (over  500  Hz),  and  geophones  (velocity 
transducers)  work  better  for  lower  frequency  waves  (less  than  500  Hz),  (Nazarian  and 
Stokoe,  1989). 


Figure  2-7.  Schematic  of  CMRP  arrangement  for  SASW  testing  (after  Nazarian  and 
Stokoe,  1989). 
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A  spectral  analyzer  and  microcomputer  are  used  to  analyze  and  record  the 
measured  data.  The  spectral  analyzer  automatically  digitizes  the  time-domain  signal,  and 
using  a  fast-Fourier  transform,  converts  the  signal  to  the  frequency  domain.  Next  the 
cross-power  spectrum  and  coherence  functions  are  generated  and  used  to  create  a  field 
dispersion  curve.  The  cross-power  spectrum  is  a  way  to  view  the  surface  wave  phase  lag 
as  a  function  of  frequency.  It  is  used  to  determine  the  surface  wave  velocity  and 
wavelength.  The  coherence  function  simply  determines  the  range  of  frequencies  over 
which  the  analyzed  signal  is  coherent,  i.e.  reasonable  and  accurate  (Nazarian  and  Stokoe, 
1989). 

Automated  SASW  testing  devices  are  still  in  the  developmental  stage,  and  as 
such,  single-unit  commercial  testing  devices  are  not  available.  Without  automation,  a 
team  of  technicians  requires  20  to  40  minutes  to  run  the  complete  series  of  tests  at  one 
site.  Thus,  the  method  is  not  considered  economically  feasible  for  general  applications 
until  it  can  be  automated. 

2.4.4  Application 

The  SASW  method  is  not  as  fast  as  the  previously  described  nondestructive  test 
methods,  but  it  does  offer  three  distinct  advantages  over  those  test  methods. 

First,  layer  thicknesses  can  be  determined  directly  from  the  testing.  With 
deflection-based  testing  and  analysis  (except  for  a  few  advanced  backcalculation 
techniques),  layer  thicknesses  are  assumed  based  on  construction  drawings,  or  core 
samples  are  drilled  to  determine  precise  thicknesses. 

Second,  the  variation  of  moduli  within  the  subgrade  can  be  determined.  Other 
tests  treat  the  entire  subgrade  as  a  single  layer  with  constant  modulus.  This  plays  a 
particularly  important  role  when  there  is  a  rigid  layer  (bedrock)  near  the  pavement 
surface.  Also,  deflection-based  testing  must  assume  the  depth  to  the  rigid  layer,  whereas 
SASW  testing  can  determine  the  depth  and  the  modulus  of  the  rigid  layer  (Nazarian  and 
Stokoe,  1989). 
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Third,  the  modulus  of  thin  layers  at  the  pavement  surface  can  be  accurately 
determined.  Large  variations  in  thin  layer  stiffness  may  only  slightly  change  the  shape  of 
a  deflection  basin,  but  the  S  ASW  method  can  easily  locate  and  determine  the  modulus  of 
thin  layers. 

2.5  Rolling  Wheel  Load 

The  newest  type  of  device  developed  for  pavement  testing  loads  the  pavement 
with  a  rolling  wheel  and  uses  sensors  to  continuously  measure  deflections.  Though 
several  loading  and  measuring  schemes  exist,  one  common  advantage  for  all  the  devices 
is  that  a  continuous  deflection  profile  is  produced.  Thus,  small  weak  zones  can  be  easily 
identified.  This  section  describes  the  operation  and  testing  of  two  different  rolling 
deflectometers. 

2.5.1  Testing  Equipment 

Rolling  weight  deflectometers  are  still  in  the  developmental  stage  and  are  not  yet 
commercially  available.  However,  some  prototypes  have  been  constructed.  Two  models 
are  described  here.  Each  loads  the  pavement  and  measures  the  deflections  using  very 
different  methods. 

2. 5. 1.1  SNRA  Rolling  Weight  Deflect ometer 

The  Swedish  Nationeil  Road  Administration  (SNRA)  began  developing  a  rolling 
weight  deflectometer  (RWD)  in  the  1980’s,  and  built  its  first  prototype  in  1991  (Amberg 
et  al.,  1991).  The  testing  device  consists  of  a  heavy  truck  with  two  transverse  laser 
profilometers.  The  profilometers  were  originally  developed  for  measuring  surface 
roughness,  but  in  this  case  were  used  to  measure  deflections. 

As  illustrated  in  Figure  2-8,  the  first  profilometer  is  mounted  a  few  feet  behind  the 
lightly  loaded  front  axle.  It  measures  the  undeflected  shape  of  the  pavement  surface.  The 
second  profilometer  is  mounted  immediately  behind  the  heavily-loaded  rear  axle  and 
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measures  the  deflected  shape  of  the  pavement.  The  difference  between  the  two  profiles 
represents  the  deflections  caused  by  the  loaded  axle. 


Figure  2-8.  Schematic  representation  of  the  SNRA  Rolling  Weight  Deflectometer. 


Preliminary  tests  with  the  prototype  indicated  that  the  measured  deflections  were 
very  repeatable  over  the  same  pavement  segment.  Tests  conducted  at  10,  30  and  50 
km/hr  (6, 19,  and  31  mph)  showed  that  deflections  markedly  decreased  as  the  vehicle 
speed  increased.  Comparisons  between  the  RWD  and  the  FWD  maximum  deflections 
(directly  behind  the  wheel  and  under  the  load  plate)  showed  reasonable  agreement  for 
wheel  velocities  of  5  km/hr.  Errors  were  attributed  to  the  differences  between  the 
dynamics  of  the  two  tests. 

2. 5. 1.2  UT  Rolling  Dynamic  Deflectometer 

Bay  et  al.  (1995)  describe  a  rolling  dynamic  deflectometer  (RDD)  constructed  at 
the  University  of  Texas  at  Austin  that  applies  a  dynamic  load  to  the  pavement  through 
two  wheels  mounted  midway  between  the  front  and  read  axles  (Figure  2-9).  The  load  is 
applied  by  a  servo-hydraulic  vibrator  and  can  generate  dynamic  vertical  forces  as  large  as 
3 10  kN  (70,000  lb).  The  servo-hydraulic  vibrator  can  generate  a  variety  of  dynamic 
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loading  functions,  including  transient  (like  the  FWD)  and  steady-state  (ranging  from  5 
to  100  Hz).  The  test  can  be  conducted  at  speeds  between  4  and  6  km/hr  (2  and  4  mph). 

The  resulting  deflections  are  measured  by  an  isolated  pair  of  wheels  resting  on  the 
pavement  directly  between  the  loading  wheels.  An  aecelerometer  records  the  vertical 
motion  of  the  wheels  to  generate  a  continuous  motion  profile.  After  some  data 
processing,  a  continuous  record  of  pavement  flexibility  can  be  constructed. 

Engine  and  Pump  for 
Servo-Hydraulic 


Receiver  Wheel  Wheel  Loading  Wheels 

Figure  2-9.  Schematic  representation  of  the  UT  Rolling  Dynamic  Deflectometer  (after  Bay  et  al., 
1995). 

Tests  have  been  completed  with  the  RDD  on  eight  flexible  pavement  sections. 
The  results  showed  that  the  test  is  particularly  good  for  identifying  weak  segments  within 
a  generally  sound  pavement  section.  Comparisons  between  the  RDD  and  FWD  found 
that  the  tests  show  similar  deflection  to  load  trends,  but  that  more  work  is  necessary  to 
fully  understand  the  relationship  between  the  two  tests  (Bay  et  al.,  1 995). 
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2.5.2  Application 

Rolling  deflectometers  show  real  promise  for  future  use  in  pavement  analysis. 
Three  aspects  stand  out  as  significant  advantages  of  rolling  load  testing.  First,  the  tests 
are  conducted  while  in  motion,  so  traffic  control  around  the  equipment  is  easier,  and  user 
(highway  travelers)  costs  are  lower.  In  the  case  of  the  SNRA  device,  high  speed  testing 
(up  to  90  km/hr  (56  mph))  may  make  the  tests  completely  transparent  to  the  user. 

Second,  the  testing  generates  a  continuous  pavement  profile.  This  helps  to  identify  weak 
segments  of  the  road  that  may  be  missed  by  conventional  stationary  tests.  Third,  for  the 
SNRA  Rolling  Weight  Defleetometer  the  applied  load  exactly  matches  traffic  loading. 
This  solves  a  problem  that  has  plagued  all  of  the  other  testing  devices. 

However,  there  are  still  several  significant  problems  to  overcome  before  rolling 
deflectometers  will  be  fully  operational.  First,  the  measurement  system  must  be 
improved  to  ensure  accurate  measurements  at  reasonable  speeds.  The  RDD  sensor  wheel 
will  likely  not  perform  well  at  even  moderate  speeds,  and  the  laser  sensors  used  by  the 
RWD  can  be  adversely  effected  by  moisture  on  the  pavement  or  a  coarse,  open-grained 
surface.  Also,  accurately  defining  the  undeflected  datum  is  difficult  because  of  small 
movements  in  the  beam  to  which  the  sensors  are  attached.  A  precisely  defined  datum  is 
essential  because  the  amplitude  of  the  measured  deflections  is  very  small  (<  1  mm). 

Second,  dynamic  interpretation  models  for  rolling  load  tests  are  not  well- 
developed.  As  previously  described,  the  dynamic  aspects  of  FWD  tests  are  frequently 
ignored  so  that  backcalculation  can  be  accomplished  more  quickly.  Rolling  loads  add  an 
extra  level  of  complexity  to  the  stationary  dynamic  load  model,  and  cause  computation 
times  to  increase  even  more.  At  this  time,  models  that  accurately  predict  the  deflection 
basin  caused  by  a  rolling  wheel  on  a  layered  pavement  system  do  not  exist. 


Chapter  3 

FLEXIBLE  PAVEMENT  DEFLECTION  MODELS 


Pavement  deflection  models  attempt  to  answer  the  question:  What  is  the  deflection 
(or  stress  or  strain)  at  location  A,  when  a  load  is  applied  at  location  B?  Many 
possibilities  exist  for  both  the  locations  and  types  of  loads  that  may  be  of  interest.  For 
example,  a  load  could  be  located  at  the  surface  or  within  a  pavement  layer,  and  the 
response  of  interest  could  be  horizontal,  vertical,  or  rotational  deflections.  Also,  applied 
loads  could  be  static  or  dynamic  over  time,  horizontal,  vertical,  or  torsional  in  direction, 
and  applied  to  a  disk,  ring,  or  square  plate. 

Ideally,  the  deflection  model  would  be  developed  to  match  the  loading  conditions 
exactly.  Inevitably,  certain  assumptions  must  be  made  to  simplify  the  model  and 
facilitate  easier  (or  more  rapid)  calculations.  For  example,  as  mentioned  in  Chapter  2, 
static  deflection  models  are  routinely  used  to  interpret  data  measured  by  dynamic  tests. 
Unfavorable  consequences  may  result  from  this  type  of  oversimplification.  This  chapter 
describes  some  flexible  pavement  deflection  models  that  have  been  developed,  and 
emphasizes  the  inherent  assumptions  of  these  models. 

3.1  Static  Models 

Static  deflection  models  predict  the  response  resulting  from  a  static  load  applied  to 
a  known  area  on  a  given  medium.  The  level  of  complexity  of  these  models  is  determined 
by  the  number  of  layers  in  the  loaded  medium,  the  stress-strain  characteristics  of  the 
materials  in  the  medium  (linear  elastic  or  non-linear  elastic),  and  the  type  of  interface 
conditions  between  individual  layers  (frictionless  or  fully-frictional). 

In  1885,  Boussinesq  developed  the  first  solutions  for  deflections  due  to  a  point  load 
applied  to  the  surface  of  a  homogeneous  linear  elastic  halfspace.  His  solutions  were 
extended  to  include  deflections  due  to  loaded  areas,  such  as  disks,  by  integrating  over  the 
loaded  area. 

The  next  advancement  in  static  deflections  occurred  in  1 943,  when  Burmister 
published  his  solution  for  displacements  in  a  two-layered  system  (a  finite  layer  over  a 
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homogeneous  halfspace).  This  allowed  pavement  analysts  to  model  the  asphalt  layer 
separately  from  the  underlying  base  and  subgrade.  Burmister  later  refined  his  solution  to 
allow  for  a  frictionless  (1945a)  or  flilly-frictional  (1945b)  interface  between  layers,  and 
for  a  three-layered  system  (1945c).  At  that  point,  pavement  analysts  could  finally  separate 
a  pavement  system  into  its  three  primary  layers:  asphalt  concrete,  base  course,  and 
subgrade. 

Linear  elastic  static  loading  models,  such  as  Burmister’ s,  make  the  following 
assumptions  about  layered  systems  (after  Huang,  1993).  A  generalized  pavement  profile 
is  illustrated  in  Figure  3-1. 

1 .  Each  layer  is  homogeneous,  isotropic,  and  linearly  elastic  with  modulus,  E, 
and  Poisson’s  ratio,  v. 

2.  The  material  is  weightless  and  infinite  in  areal  extent. 

3.  Each  layer  has  a  finite  thickness,  h,  but  the  lowest  layer  is  infinite  in 
thickness. 

4.  A  uniform  pressure,  q,  is  applied  on  the  surface  of  the  pavement  over  a 
circular  area  of  radius  a. 

5.  Continuity  equations  are  satisfied  at  layer  interfaces,  as  indicated  by  the  same 
vertical  stress,  shear  stress,  vertical  displacement,  and  radial  displacement. 


r:Tra'' 

Layer  1  E,,  v, 

\ 

'hi 

/ 

Layer  2  Ez,  V2 

/ 

< 

L _ _ _ _ _ _ _ _ V. 

v!" . 

’  ^  7 

Layer  n  E„,  v„ 

\ 

00 

'if 


Figure  3-1.  Illustration  of  an  n-layer  system  subjected  to  a  disk  load  (after  Huang, 
1993). 
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The  first  assumption,  which  states  that  all  materials  are  linear  elastic,  is  a  very  important 
feature  of  these  models.  The  linear  aspect  allows  the  use  of  superposition  to  calculate 
deflections  caused  by  multiple  loads.  For  example,  the  deflection  at  the  midpoint 
between  two  equally  loaded  tires  is  equal  to  twice  the  magnitude  of  deflection  at  the  same 
distance  from  a  single  tire.  The  elastic  portion  of  the  assumption  recognizes  that  moving 
traffic  loads  cause  small  deformations  that  are  recoverable,  even  though  subgrade  soils  do 
not  behave  elastically  when  subjected  to  large,  stationary  loads  (Huang,  1993). 

These  assumptions  can  be  applied  to  a  multi-layered  system,  and  with  the  use  of 
computers,  programs  have  been  developed  that  quickly  solve  for  stresses,  strains,  and 
displacements  on  the  surface  and  within  the  pavement  system  (Huang,  1993).  ELSYM5, 
KENLAYER,  WESLEA,  BISAR,  and  EVERSTRESS  all  use  layered  linear  elastic  theory 
to  solve  the  static  load  problem. 

Non-linear  elastic  models  recognize  that  pavement  system  materials,  especially 
subgrade  soils,  do  not  have  a  constant  modulus  as  stresses  and  strains  increase.  Non¬ 
linear  models  require  knowledge  of  a  soil’s  constitutive  relationship  (Huang,  1993), 
which,  in  the  simplest  case,  relates  the  stiffness  of  a  soil  element  to  its  state  of  stress. 
When  a  pavement  system  is  modeled  as  a  combination  of  many  small  elements,  a 
constitutive  relationship  allows  the  moduli  to  be  adjusted  for  each  element  within  a  soil 
layer,  depending  on  the  overburden  and  applied  stress. 

Chang  et  al.  (1992)  studied  the  potential  for  non-linear  behavior  in  the  FWD  test. 

By  varying  the  applied  loads  at  a  single  site,  the  deflection  behavior  was  obtained.  Chang 
et  al.  determined  that  non-linearities  were  especially  pronounced  where  the  subgrade  was 
soft  or  the  asphalt  pavement  was  relatively  thin.  In  these  cases,  non-linear  effects 
increased  peak  deflections  by  as  much  as  50%  under  very  heavy  FWD  loads  (88.96  kN 
(20,000  lb)).  However,  deflections  more  than  3  ft  away  from  the  load  plate  were  nearly 
unchanged. 

Several  computer  programs  have  been  developed  to  model  the  non-linear  behavior 
of  multi-layered  systems.  MICH-PAVE,  developed  at  Michigan  State  University,  and 
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ILLI-PAVE,  developed  at  the  University  of  Illinois,  both  use  finite  element  methods  to 
calculate  stresses,  strains,  and  deflections  due  to  static  loads  on  non-linear  pavement 
systems. 

3.2  Dynamic  Models 

Many  pavement  analysts  use  static  models  to  interpret  the  data  obtained  from 
dynamic  tests,  such  as  the  FWD  or  Road  Rater.  Although  this  practice  is  common, 
certain  conditions  can  result  in  gross  misinterpretation  of  the  data.  For  example,  Roesset 
and  Shao  (1985)  concluded  that  a  static  interpretation  of  Dynaflect  (steady-state  dynamic 
load)  deflections,  when  the  soil  is  underlain  by  a  stiff  layer  (bedrock)  at  depths  of  less 
than  60  feet,  can  result  in  substantial  errors  (up  to  70%)  in  pavement  layer  moduli.  The 
errors  were  attributed  to  dynamic  amplifications  at  certain  loading  frequencies.  Dynamic 
effects  are  not  as  pronounced  in  FWD  tests  because  the  impact  load  excites  a  broad  range 
of  frequencies.  However,  the  AC  moduli  still  could  be  underpredicted  by  as  much  as 
40%  when  the  FWD  dynamic  response  is  not  considered. 

This  section  describes  dynamic  deflection  models  that  use  a  variety  of  mathematical 
techniques.  First,  models  that  consider  a  dynamic  load  applied  at  one  location  are 
explained.  Next,  moving  load  models  are  described. 

3.2.1  Stationary  Load  Models 

In  this  thesis,  dynamic  stationary  load  models  refer  to  those  in  which  the  load  is 
applied  at  the  same  location  on  a  pavement  throughout  the  duration  of  loading.  Dynamic 
steady-state  and  FWD  tests  both  fall  into  this  category. 

3. 2. 1.1  Linear-Elastic  Dynamic  Models 

Determining  the  response  of  a  soil  deposit  to  dynamic  loading,  whether  caused  by 
explosions,  moving  traffic,  or  earthquakes,  falls  into  the  mathematical  realm  of  wave 
propagation.  In  1904,  Lamb  first  published  the  solution  for  waves  passing  through  an 
elastic  halfspace  (Meier,  1995).  Thomson  (1950)  was  the  first  to  develop  the  solution  for 
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waves  passing  through  a  layeredhsXfs^QCQ.  Thomson’s  solution  used  matrices  that 
linearly  related  stresses  at  the  top  of  a  layer  to  the  stresses  at  the  bottom,  and  ensured  that 
particle  velocities  at  the  bottom  of  one  layer  were  identical  to  those  at  the  top  of  the 
underlying  layer.  Using  Thomson’s  method,  the  stresses  at  the  top  or  bottom  of  any  layer 
could  be  related  to  the  stresses  at  the  surface. 

Haskell  (1953)  used  Thomson’s  matrix  formulation  to  generate  the  phase  velocity 
dispersion  equations  for  elastic  surface  waves  (Rayleigh  and  Love  waves)  through  a 
multi-layered  halfspace.  Haskell  also  corrected  an  error  in  Thomson’s  derivation,  thus 
his  results  produced  what  is  now  called  the  Haskell-Thomson  transfer  matrix  method 
(Meier,  1995). 

Kausel  and  Roesset  (1981)  used  the  Haskell-Thomson  transfer  matrix  approach  to 
derive  layer  stiffiiess  matrices  that  relate  applied  loads  to  displacements.  The  stiffhess 
matrix  solution  is  not  more  general  than  the  Haskell-Thomson  transfer  matrix  approach, 
but  it  is  more  efficient  for  numerical  integration.  When  layer  thicknesses  are  small 
compared  to  the  wavelengths  of  interest,  the  transcendental  functions  which  govern 
vertical  displacements  can  be  linearized.  Algebraic  functions  then  replace  the 
transcendental  functions,  and  the  eigenvalue  problems  for  wave  propagation  can  be 
solved  by  standard  techniques.  This  results  in  significant  gains  in  computational 
efficiency. 

Using  this  stiffness  matrix  approach,  Kausel  and  Peek  (1982)  developed  an  explicit, 
closed-form  solution  for  functions  (Green’s  fimctions  -  described  in  Chapter  4)  that  relate 
displacements  in  a  layered  strata  to  unit  harmonic  loads.  Solutions  were  developed  for 
several  loading  schemes  including  vertical  disk  loads,  horizontal  point  loads,  torsional 
ring  loads,  and  other  permutations  of  load  direction  and  areas.  Kausel  created  a  computer 
program  called  PUNCH,  that  uses  this  method  of  calculating  the  Green  functions. 
PUNCH  has  been  used  by  several  researchers  to  develop  dynamic  pavement  deflection 
models. 

Davies  and  Mamlouk  (1985)  studied  deflection  basins  generated  by  dynamic 
steady-state  testing  devices,  and  concluded  that  the  basins  were  substantially  different 
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than  those  produced  by  static  loads.  The  application  of  static  deflection  analyses 
applied  to  dynamically  produced  basins  may  yield  misleading  results  if  the  operating 
frequency  is  close  to  the  resonant  frequency  of  the  pavement  system.  Other  researches 
have  studied  the  influence  of  a  rigid  layer  (bedrock)  at  a  shallow  depth  (less  than  60  ft)  on 
deflection  basins.  Roesset  and  Shao  (1985)  showed  that  static  analysis  of  the  basin 
produced  by  steady-state  (Dynaflect)  and  transient  (FWD)  dynamic  loads  can  lead  to 
large  errors  in  backcalculated  layer  moduli. 

Sanchez-Salinero  (1987)  used  PUNCH  to  model  waves  produced  during  seismic 
testing  of  soils.  After  conducting  a  series  of  parametric  studies,  he  developed  rules  for 
discretizing  a  layered  system  to  yield  the  most  accurate  solution. 

Ketchum  (1993)  developed  a  dynamic  model  that  uses  PUNCH  to  interpret  FWD 
tests  on  airport  pavements.  He  focused  on  the  effects  of  a  frozen  layer  at  shallow  depth 
on  the  dynamic  response  of  pavements,  and  weis  able  to  formulate  a  backcalculation 
procedure  to  determine  layer  moduli.  Ketchum  found  that  pavements  could  exhibit  a 
significant  resonant  amplification  at  certain  frequencies  when  a  rigid  layer  is  at  a  shallow 
depth. 

To  generate  dynamic  deflection  basins  for  training  a  neural  network  backcalculation 
scheme,  Meier  (1995)  used  Sanchez-Salinero’s  version  of  PUNCH.  He  optimized  the 
code  to  calculate  only  vertical  deflections  due  to  a  vertical  disk  load,  and  added  a 
subroutine  that  automatically  discretized  the  pavement  layers  into  properly-sized 
sublayers. 

3. 2. 1.2  Viscoelastic  Dynamic  Models 

Asphalt  concrete  is  considered  a  viscoelastic  material.  That  is,  the  stress  in  the 
eisphalt  is  a  function  of  the  strain  and  the  strain  rate.  Viscoelastic  materials  are  often 
represented  in  models  by  a  linked  combination  of  a  linear  spring  and  dashpot  —  the  stress 
in  the  spring  is  directly  proportional  to  the  strain,  whereas  the  stress  in  the  dashpot  is 
proportional  to  the  strain  rate.  For  example,  the  Maxwell  model,  illustrated  in  Figure  3-2, 
is  represented  by  a  spring  and  dashpot  in  series.  For  a  constant  applied  stress,  an 


instantaneous  elastic  extension  of  the  spring  occurs.  At  the  same  time,  the  dashpot 
begins  to  extend,  and  over  a  period  of  time  the  strain  will  steadily  increase  (Haddad, 
1995). 
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Figure  3-2.  Illustration  of  the  properties  of  a  viscoelastic  Maxwell  model  (after  Haddad, 
1995). 


Because  FWD  and  moving  wheel  loads  are  inherently  time-dependent  loads,  it  is 
logical  to  build  a  viscoelastic  model  to  predict  deflections.  For  example,  KENLAYER 
has  the  option  to  consider  dynamic  loads  by  assigning  viscoelastic  properties  to  each 
layer  (Huang,  1993).  SAPSI,  developed  at  the  University  of  California  at  Berkeley,  uses 
the  finite  element  method  with  viscoelastic  properties  to  calculate  deflection  responses 
resulting  from  multiple  dynamic  surface  loads  (Chatti  et  al.,  1995). 

3.2.2  Rolling  Load  Models 


Relatively  few  models  have  been  created  to  predict  the  displacements  from  a 
moving  load  on  a  pavement  system.  Huang  (1993)  describes  viscoelastic  models  using 
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KENLAYER  and  VESYS  to  predict  the  time  history  of  deflection  at  a  single  location 
as  the  wheel  passes  over  it.  For  these  models,  a  series  of  individual  loads  are  applied  to 
the  point  at  each  time  increment.  The  amplitude  of  each  point  load  increases  and 
decreases  following  the  pattern  of  a  haversine.  The  duration  of  the  combined  pulse  is 
dependent  on  the  speed  of  the  vehicle  and  the  tire  contact  area. 

An  alternative  approach  that  is  particularly  applicable  to  railroads  is  to  treat  the 
asphalt  surface  layer  as  a  beam  fully  supported  by  an  elastic  or  viscoelastic  foundation. 
Achenbach  and  Sun  (1965)  used  a  viscoelastic  model  to  study  the  dynamic  response  of  a 
beam  resulting  from  a  rapidly  moving  load.  Holder  and  Michalopoulos  (1977) 
investigated  the  response  of  a  beam  supported  by  a  Winkler-type  elastic  foundation  with 
inertia  (i.e.  not  massless)  and  damping  subjected  to  a  moving  load.  Steady-state  solutions 
for  beam  deflections  were  developed;  that  is,  the  deflections  an  observer  moving  with  the 
load  would  see. 

Trochanis  et  al.  (1987)  predicted  the  beam-foundation  contact  pressure,  beam 
deflection,  and  bending  moment  resulting  from  a  moving,  oscillating  load  (train  wheel 
set)  on  a  fully-supported  beam.  The  system  was  modeled  as  an  infinitely  long  beam 
supported  by  an  arbitrary  linear-elastic,  damped  foundation.  The  model  showed  that  for  a 
moving  load,  the  maximum  deflection  occurs  slightly  behind  the  wheel,  as  illustrated  in 
Figure  3-3. 

Deflection 


Figure  3-3.  Static  and  dynamic  deflections  resulting  from  a  rolling  wheel  load  (Trochanis  et  al.,  1987). 


Chapter  4 

DYNAMIC  LAYERED  ELASTIC  ANALYSIS  USING  GREEN’S  FUNCTIONS 


The  previous  chapter  described  a  number  of  different  types  of  pavement 
deflection  models.  This  chapter  will  focus  on  one  particular  type  --  a  dynamic  layered 
elastic  model.  In  this  research,  a  dynamic  layered  elastic  model  was  developed  using  the 
Green’s  function  approach.  This  chapter  describes  the  Green’s  function  concept, 
applications,  and  the  specific  solutions  developed  for  this  project. 

4.1  Green’s  Functions 

A  Green’s  function  is  a  mathematical  tool  that  can  be  used  to  find  solutions  to 
many  boundary  value  problems  associated  with  either  ordinary  or  partial  differential 
equations.  The  Green’s  function  concept  was  first  introduced  by  G.  Green  in  1828,  but 
its  use  as  a  practical  analysis  tool  did  not  begin  to  develop  until  the  middle  of  this  century 
(Roach,  1970).  The  Green’s  function  itself  is  unique  for  each  problem  and  set  of 
boundary  conditions.  However,  once  the  Green’s  function  is  determined,  it  reduces  the 
solution  to  an  integral  equation  with  boundary  conditions,  rather  than  a  differential 
equation.  For  numerical  analysis,  the  integral  form  is  easier  to  solve  than  the  differential 
equation  (Roach,  1970). 

The  Green’s  function  approach  has  been  applied  to  many  practical  applications  in 
engineering.  For  example.  Roach  (1970)  used  Green’s  functions  to  solve  boundary  value 
problems  associated  with  the  wave  equation  and  the  diffusion  equation.  Trim  (1990) 
described  the  use  of  Green’s  functions  to  find  the  displacement  in  a  horizontal,  taut  string 
that  is  fixed  at  one  end  while  the  other  end  is  moved  harmonically  in  the  vertical 
direction. 

Hartzell  (1978)  described  a  Green’s  function  approach  for  predicting  earthquake 
grovmd  motion.  An  earthquake  was  modeled  eis  a  series  of  point  ruptures  on  a  fault  zone. 
A  total  earthquake  motion  was  produced  by  superimposing  the  phase-delayed  ground 
motion  resulting  from  each  point  rupture.  Hartzell  scaled  previously  recorded  aftershocks 
(small  tremors)  to  serve  as  the  response  at  each  point  source  —  these  small  motions  are 
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empirical  Green’s  functions.  The  Green’s  function  method  of  dividing  a  large  problem 
into  the  sum  of  smaller,  easier  to  solve  problems  is  a  common  feature  of  Green’s  function 
analyses,  and  was  employed  in  this  thesis. 

4.2  Dynamic  Layered  Elastic  Model  using  Green’s  Functions 

In  this  research,  the  Green’s  function  can  be  viewed  as  a  simple  transfer  function, 
GfcOj,  R),  that  relates  displacements,  f/fco;,  R),  to  load, 

U((o,,R)=G(co,,R)*F(co,) 

Combining  the  displacements  at  all  frequencies  produces  the  desired  displacement  time 
history. 

+00 

i=—co 

This  formulation  was  used  to  build  a  pavement  deflection  model,  which  was  coded  into 
the  computer  program  BASIN. 

The  following  sections  describe  the  manner  in  which  the  Green’s  functions 
approach  was  applied  in  BASIN.  First,  conversion  from  the  time  domain  to  the 
frequency  domain  for  the  analysis  of  transient  loads  and  responses  is  described.  Next,  the 
generalized  method  for  calculating  the  elastodynamic  Green’s  functions  is  explained. 
Finally,  the  last  sections  illustrate  the  schemes  used  to  represent  both  FWD  and  rolling 
wheel  loads. 

4.2.1  Frequency  Domain  Analysis 

As  mentioned  in  Chapter  3,  the  Green’s  functions  solved  by  Kausel  and  Peek 
(1982)  describe  the  harmonic  deflections  produced  by  steady-state  harmonic  unit  loads. 
To  use  these  Green’s  functions  for  calculating  deflections  resulting  from  transient  loads, 
the  transient  load  must  be  represented  in  the  frequency  domain  as  a  combination  of 
harmonic  loads  of  different  amplitude,  phase,  and  frequency.  This  transformation  and 
combination  approach  is  commonly  known  as  Fourier  superposition  analysis. 
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Figure  4-1  illustrates  a  general  schematic  for  analysis  in  the  frequeney  domain. 

In  the  first  step,  the  load  time  history  is  discretized  into  individual  points  spaced  at  a 
constant  time  interval.  The  discretized  loading  is  then  decomposed  into  the  sum  of  a 
series  of  simple  harmonic  loads  using  a  discrete  Fourier  transform.  The  Green’s 
functions  are  then  solved  to  determine  the  amplitude  and  phase  of  harmonic  displacement 
resulting  from  a  unit  load  at  each  frequency.  Next,  the  Green’s  functions  are  multiplied 
by  the  Fourier  series  of  the  loading  to  determine  the  Fourier  series  of  the  pavement 
defleetion.  Finally,  the  sum  of  the  harmonic  responses  are  transformed  back  into  the  time 
domain  by  an  inverse  discrete  Fourier  transform.  This  method  of  superimposing  the 
harmonic  responses  is  only  valid  for  linear-elastic  materials. 

4.2.2  Green’s  Function  Solution 

As  described  previously,  an  explicit,  closed-form  solution  for  the  Green’s 
functions  corresponding  to  dynamic  loads  acting  on  layered  strata  was  used  to  build  the 
pavement  deflection  model  BASIN.  This  seetion  very  briefly  summarizes  the  general 
solution  process.  The  complete  derivation  can  be  found  in  Kausel  and  Peek  (1982). 
Additionally,  an  excellent  overview  is  given  by  Meier  (1995). 

The  solution  developed  by  Kausel  and  Peek  (1982)  uses  the  following  steps  to 
solve  for  harmonic  displacements  in  a  layered  strata.  First,  the  stiffness  matrix  for  each 
layer  must  be  computed.  Each  stiffness  matrix  is  a  function  of  the  material  properties  and 
thickness  of  the  layer.  Additionally,  the  equation  for  the  stiffness  matrix  is  dependent  on 
the  loading  frequency  and  wave  number.  The  wave  number  is  obtained  by  a  spatial 
Fourier  transformation  of  the  loading  function  (this  is  analogous  to  a  time-dependent 
ftinction  transformed  into  the  frequency  domain)  and  is  used  to  include  the  distance 
between  the  loaded  area  and  the  observation  point  as  part  of  the  solution.  A  global 
stiffhess  matrix,  K,  is  obtained  by  overlapping  the  individual  layer  matrices  at  their  layer 
interfaees.  The  global  stiffiiess  matrix  relates  the  displacements,  U  ,to  the  loads,  P ,  i.e. 


P  =  KU 


Response  represented  as  sum  of  a  Summation  of  harmonic  responses 

series  of  harmonic  responses  to  produce  time  history  response 

Figure  4-1.  Schematic  illustration  of  Fourier  superposition  analysis  using  Green’s  functions  (after  Kramer, 
1996). 

By  formally  inverting  the  stiffness  matrix,  an  expression  that  solves  directly  for  the 
displacements  is  foimd, 

U  =  K'P  =  FP 
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where,  F  =  K'^  represents  the  flexibility  matrix.  After  the  flexibility  matrix  is 
established,  it  is  multiplied  by  the  loading  vector  to  calculate  the  displacements. 

Kausel  and  Peek  (1982)  developed  solutions  for  the  displacements  produced  by 
several  types  of  loading.  Displacements,  in  cylindrical  coordinates,  are  calculated  in  the 
vertical,  radial,  and  tangential  directions.  Loads  can  be  applied  horizontally  or  vertically 
to  a  disk,  ring  or  point;  additionally,  torsional  and  rocking  loads  are  solved. 

4.2.3  FWD  Load  Representation 

The  load  pulse  produced  by  an  FWD  was  developed  to  approximate  the  deflection 
pulse  produced  by  a  moving  truck  load  (Hoffman  and  Thomson,  1982).  Barksdale  (1971) 
reported  the  actual  load  pulse  caused  by  a  moving  wheel  closely  resembled  a  haversine, 
illustrated  in  Figure  4-2.  Example  pulses  from  several  FWD  devices  are  shown  in  Figure 
4-3.  The  shape  and  duration  of  the  pulse  are  controlled  by  the  mass,  height  of  drop,  and 
buffer  system  stif&iess. 


Figure  4-2.  Typical  wheel  load  applied  by  a  moving  truck 
(Barksdale,  1971). 


For  dynamic  analysis  of  pavement  systems,  it  is  computationally  easier  to  model 
the  load  as  a  smooth  function.  The  level  of  agreement  between  these  functions  and  the 
“real”  load  shapes  are  most  easily  distinguished  when  the  loads  are  represented  in  the 
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frequency  domain.  Figure  4-4  shows  time-domain  loading  and  the  Fourier  amplitude 
spectrum  for  a  triangle,  haversine  and  real  pulse,  all  with  equal  peak  amplitude  and 
duration.  Because  of  its  sharp  peak  and  straight  sides,  the  triangular  pulse  has  more 
energy  in  the  100  to  200  Hz  range  than  the  haversine  or  measured  pulse  (Meier,  1995). 
Huang  (1993)  and  Ong  et  al.  (1991)  suggested  that  the  ideal  approximating  shape  is  a 
haversine  —  that  shape  was  used  for  all  analyses  in  this  project.  There  is  an  additional 
benefit  to  using  a  haversine  rather  than  a  triangular  pulse  in  this  project.  Because  the 
haversine  function  has  less  high  frequency  energy  than  the  triangular  wave,  a  smaller 
number  of  discrete  frequencies  are  necessary  to  accurately  define  the  load.  The  upper 
boimd  of  the  range  of  frequencies  used  to  characterize  the  loading  is  called  the  cutoff 
frequency  —  i.e.  energy  at  frequencies  higher  than  the  cutoff  is  neglected.  The  range  of 
frequencies  between  0  Hz  (static)  and  the  cutoff  frequency  is  also  known  as  the 
bandwidth.  Because  GREEN  calculates  a  solution  for  each  frequency  of  interest,  a  lower 
cutoff  frequency  requires  less  computation  time  to  determine  the  response  to  the 
haversine  loading. 

Meier  (1995)  showed  that  a  haversine  constructed  to  match  impulse  and  duration 
of  a  measured  pulse  also  matched  over  the  range  of  0-75  Hz  in  the  frequency  domain 
(Figure  4-5).  Differences  are  apparent  at  frequencies  above  75  Hz,  but  the  amplitudes  are 
very  small  relative  to  those  at  lower  frequencies. 

Because  computational  times  for  GREEN  increase  in  direct  proportion  to  the 
number  of  frequencies  needed  to  represent  the  loading  pulse,  it  is  desirable  to  minimize 
the  cutoff  frequency  and  maximize  the  increment  between  frequencies.  To  examine  the 
effects  of  limiting  the  cutoff  frequency,  a  FFT  was  performed  on  a  1024-point  time 
reeord  with  a  30  msec  (60  point)  haversine  pulse  at  the  beginning.  The  30  msec  duration 
was  selected  as  t5q)ical  for  FWD  tests.  The  total  frequency  range  for  this  record  was  0  to 
996.45  Hz,  with  511  discrete  frequencies  spaced  at  1.95  Hz.  To  limit  computation  time 
the  cutoff  frequency  was  set  at  76.05  Hz,  therefore  only  the  first  40  diserete  frequencies 
were  used.  An  inverse  FFT  on  the  cutoff  (or  bandwidth-limited)  spectrum  was  then  used 
to  eonvert  the  record  back  to  the  time  domain.  Figure  4-6  shows  that  there  is  virtually  no 
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difference  between  the  original  and  bandwidth-limited  pulses.  Meier  (1995)  produced 
similar  results  using  frequencies  below  75.47  Hz  with  a  2.36  Hz  interval.  Section  7.2.2 
describes  the  sensitivity  of  the  deflection  basin  to  the  cutoff  frequency  and  shows  that  this 
bandwidth  is  acceptable  for  accurately  calculating  deflection  basins. 


KUAB  FWD  DYNATEST  HWD 


DYNATEST  FWD 


PHONIX  FWD 


TIME,  msec 


Figure  4-3.  Typical  load  pulses  for  four  different  FWD  devices  (Bentson  et  al.,  1989). 


Load  Factor 


4.2.4  Rolling  Wheel  Load  Representation 

To  model  the  load  applied  by  a  rolling  wheel,  it  is  necessary  to  develop  a  method 
that  can  approximate  a  moving  load  applied  over  time  and  space.  Ideally,  the  load  would 
be  represented  by  a  sequence  of  infinitesimally-spaced  point  loads.  However,  this  is 
computationally  impossible.  Instead  the  line  load  can  be  achieved  by  discretizing  the 
load  into  a  series  of  pulses,  such  that  the  combination  of  pulses  approximates  the  desired 
load  in  time  and  space. 

Four  components  of  the  individual  pulses  must  be  considered  to  develop  the 
rolling  wheel  representation:  the  area  on  which  the  load  is  applied,  the  shape  of 
individual  load  pulses,  the  duration  of  each  pulse,  and  the  spacing  between  each  pulse. 
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Figure  4-6.  Original  and  bandwidth-limited  haversine  load  pulses  (after  Meier,  1995). 

4.2.4. 1  Loaded  Area 

The  loaded  area  was  modeled  as  a  circle  with  a  radius  that  is  a  function  of  the  tire 
pressure  and  applied  load.  For  example,  a  radial  tire  inflated  to  620.6  kPa  (90  psi) 
carrying  a  load  of  30,000  N  (6,745  lb)  has  a  contact  area  of  48,387  mm^  (75  in^).  A 
circular  disk  with  this  area  would  have  a  radius  of  120  mm  (4.7  in). 

4.2.4.2  Load  Pulse  Shape 

Figure  4-7  illustrates  the  maimer  in  which  symmetrical  load  pulses  can  be  applied 
in  succession  to  approximate  a  moving  point  load  of  constant  amplitude.  In  this  figure, 
each  individual  pulse  has  a  small  amplitude  with  a  0. 1  sec  duration,  and  is  represented  by 
a  dashed  line.  The  bold  line  above  the  individual  pulses  is  calculated  by  summing  the 
load  contribution  of  each  pulse  at  each  time  increment.  From  this  figure,  it  is  clear  that  a 
simple,  symmetrical  function  is  most  desirable  for  modeling  the  load  pulse.  For  reasons 


previously  described,  a  haversine  was  selected  as  the  standard  load  pulse  shape  for 
modeling  moving  loads. 
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Figure  4-7.  Illustration  of  method  for  combining  different  load  pulse  shapes  to  achieve  a  common  line 
load. 


4.2.4. 3  Load  Pulse  Duration 

The  duration  of  the  load  pulse  from  a  moving  wheel  influences  the  nature  of  the 
pavement  deflections.  Two  alternatives  were  examined  in  the  selection  of  an  appropriate 
load  pulse  duration. 

Barksdale  (1971)  developed  a  chart  for  approximating  the  duration  of  applied  load 
for  trucks  moving  at  several  speeds,  ranging  from  1  to  45  mph,  illustrated  in  Figure  4-8. 
The  charts  were  derived  from  an  elastic  finite  element  model  using  linear  elastic 
elements.  Barksdale’s  chart  was  generalized  to  calculate  stress  pulse  times  at  all  depths 
beneath  a  pavement  surface.  However,  for  this  project,  only  the  pulse  duration  at  the 
pavement  surface  is  of  interest.  For  a  vehicle  speed  of  50  km/hr  (31  mph),  Barksdale 
predicts  a  30  msec  pulse  duration. 

Alternatively,  the  duration  could  be  decided  by  the  length  of  time  that  the  tire  is  in 
contact  with  one  point  on  the  pavement  surface  at  a  specific  speed.  Continuing  with  the 
example  begun  in  Section  4.2.4.1,  if  the  wheel  were  rolling  at  50  km/hr  =  13.89  m/s,  then 
the  time  that  the  tire  is  in  contact  with  any  given  point  can  be  calculated  as  follows. 


Diameter  _  0.24m 
Velocity  13.89%,. 


0.0 173  sec 


1 7.3msec 
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With  this  scheme  the  pulse  duration  changes  when  the  speed,  tire  pressure,  or  load 
changes. 


Q  4  B  12  16  20  24  28 


DEPTH  BENEATH  PAVEMENT  SURFACE  (INCHES) 

Figure  4-8.  Variation  of  equivalent  principal  stress  pulse  time  with  vehicle 
velocity  an(d  (depth  (Barksdale,  1971). 

For  the  prototype  model  which  only  allows  constant  load  and  velocity,  the  pulse 
durations  reported  by  Barksdale  (1971)  are  used.  Improved  algorithms  that  allow 
variable  load  and  velocity  could  conceivably  use  a  scheme  to  automatically  calculate  the 
pulse  duration  based  on  the  tire  contact  duration. 

For  the  selected  pulse  duration,  a  constant  line  load  can  be  achieved  by  adjusting 
the  amplitude  of  each  individual  pulse,  as  illustrated  by  Figure  4-9. 
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Load  Load  Load 


Figure  4-9.  Illustration  of  method  for  using  different  load  pulse  duration  and  amplitude  to  achieve  a 
common  line  load. 


4.2A.4  Load  Pulse  Spacing 

The  final  component  of  the  rolling  load  model  is  the  spacing  between  individual 
load  pulses.  The  same  line  load  can  be  achieved  with  different  spacing  between  pulses  of 
the  same  duration,  but  different  magnitude,  as  illustrated  by  Figure  4-9.  Individual  load 
pulses  are  represented  by  thin,  dashed  lines.  The  bold  line  illustrates  the  total  applied 
load,  and  is  equal  to  the  sum  of  the  pulses  at  each  time  increment.  In  each  case,  the 
amplitude  of  an  individual  pulse  is  calculated  using  the  following  formula: 


Pulse  Amplitude  =  Line  Load  Magnitude  * 


2  *  Pulse  Spacing 
Pulse  Duration 


To  improve  required  computation  time,  a  minimum  number  of  pulses  are 
desirable.  However,  the  spacing  between  pulses  has  a  significant  impact  on  the 
calculated  deflections.  Increasing  the  load  pulse  spacing  necessarily  increased  the 
distance  traveled  by  the  wheel  between  load  pulses.  Therefore,  deflection  responses  from 
more  widely-spaced  pulses  contribute  to  the  total  deflection.  As  the  load  pulse  spacing 
increases,  the  model  less  closely  approximates  a  moving  wheel  as  a  sequence  of 
infinitesimally-spaced  point  loads.  Chapter  5  more  completely  describes  this  concept. 
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and  Chapter  7  contains  the  details  of  a  parametric  study  that  determined  the  optimum 
spacing  between  load  pulses. 
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Figure  4-10.  Illustration  of  method  for  using  different  load  pulse  amplitude  and  spacing  to  achieve  a 
common  line  load. 


Chapter  5 

BASIN  CONCEPTUAL  DEVELOPMENT 


This  chapter  describes  an  approach  developed  to  compute  deflection  basins  using 
Green’s  functions.  Initially,  a  simple,  simulated  deflection  response  function  was  used  in 
place  of  the  Green’s  function  to  develop  the  logic  process  for  deflection  basin 
calculations.  After  the  basic  logic  had  been  developed  and  verified,  the  complex- valued 
Green’s  function  response  solved  by  Kausel  and  Peek  (1982)  was  used  to  compute 
deflection  basins  for  layered  elastic  systems.  This  final  program  was  coded  in  the 
program  BASIN. 

5.1  Algorithm  Development 

The  algorithm  for  computing  deflection  basins  was  developed  by  working  through 
several  loading/observation  scenarios  -  from  simple  FWD  loading  to  complex,  dual 
rolling  wheel  loading.  Initially,  a  simple  deflection  response  function  was  developed  to 
simulate  the  time  history  of  deflection  due  to  a  unit  impulse  load  at  a  given  distance  from 
an  observation  point.  This  simple  function  was  used  to  develop  the  logic  for  computing 
observation  point  deflections  caused  by  time- varying  loads  from  both  stationary  and 
moving  sources.  This  section  presents  the  general  concepts  used  to  compute  deflections 
using  a  Green’s  function  approach. 

5.1.1  Simulated  Green’s  Function 

A  simulated  Green’s  function  response  to  a  unit  impulse  load  was  developed  as  a 
function  of  time  and  distance  from  the  impact  point.  The  function  was  designed  so  that 
the  deflections  observed  at  nearby  observation  points  would  arrive  quickly  and  with  large 
amplitude.  Responses  farther  from  the  load  would  arrive  later  and  with  smaller  amplitude 
to  simulate  the  effects  of  radiation  and  material  damping.  Additionally,  a  small 
amplitude,  low  frequency,  transient  vibrating  response  was  included  in  the  simulated 
deflection  response.  This  generates  a  small  rebound  (positive  deflection)  after  the  initial 
deflection  at  small  measurement  distances  from  the  load. 
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The  function  developed  for  the  simulated  deflection,  D(R,T),  was  manipulated 
to  produce  a  plausible  response  shape  (i.e.  coefficients  were  added  and  adjusted  until  the 
shape  looked  reasonable).  The  deflection,  D,  is  a  function  of  time,  T,  and  distance,  R,  and 
is  of  the  form: 


D{R, T)  =  e  *C*(a*b  +  e  *  A sin(cor)) 


where: 
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Figure  5-1  illustrates  the  simulated  time  history  of  deflection  measured  at  several 
distances  from  the  load.  Units  have  deliberately  been  omitted  from  all  the  figures  in  this 
section  to  emphasize  that  these  are  simulated  responses. 
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The  simulated  response  function  was  used  while  developing  the  logic  required 
for  prediction  of  observation  point  deflections  under  a  variety  of  loading  conditions. 
First,  this  function  was  used  to  compute  deflection  basins  produced  by  falling  weight 
deflectometer  (FWD)  tests  —  a  case  in  which  the  load  varies  with  time  but  not  position 
(i.e.  the  load  is  stationary)  and  the  observation  point  is  stationary.  Next,  the  Green’s 
function  approach  was  applied  to  the  case  of  a  single  moving  load  of  constant  amplitude 
and  velocity  and  passing  by  a  stationary  observation  point.  Finally,  this  approach  was 
extended  to  predict  the  shape  of  a  deflection  basin  that  could  be  measured  by  deflection 
sensors  moving  with  the  load. 


Deflection 


Figure  5-1.  Simulated  deflection  response  caused  by  a  unit  impulse  load  for  several 
distances  from  the  impact  site. 

5.1.2  Case  1:  Falling  Weight  Deflectometer 

As  described  in  Chapter  2,  commonly  only  the  peak  deflections  at  each  sensor  are 
recorded  in  a  FWD  test.  Alternatively,  the  system  could  record  the  time  history  of 
deflection  at  each  sensor  to  capture  the  dynamic  effects  of  the  load. 

Case  1  models  a  FWD  test  to  determine  the  time  history  of  deflection  at  several 
sensors  located  at  different  distances  from  the  load.  Figure  5-2  shows  the  arbitrary  shape 


55 


of  the  load  used  in  this  process.  The  continuous  actual  loading  curve  was 
approximated  by  a  series  of  point  loads,  P,  separated  by  an  equal  time  step,  At.  Each  of 
these  point  loads  was  assigned  an  impulse  value,  I(i)  =P(i)  *  At,  that  was  multiplied  by  the 
simulated  deflection  response  for  that  distance.  Figure  5-3  shows  the  response  at  R(10) 
due  to  each  impulse;  the  response  due  to  the  point  load  P(i)  (from  Figure  5-2)  is 
highlighted.  To  find  the  total  response,  the  principle  of  superposition  was  used  to 
combine  the  deflections  at  each  time  step.  The  resulting  total  response  at  this  sensor  is 
shown  in  Figure  5-4. 


Figure  5-2.  Simulated  impact  load  function  conversion  to  discrete  point  loads. 


This  process  was  repeated  for  each  measurement  distance.  Figure  5-5a  shows  the 
time  history  of  deflection  at  several  distances.  Using  the  peak  values  from  each  of  the 
curves,  the  deflection  basin  measured  by  a  typical  FWD  was  created  (Figure  5-5b)  Note 
that  the  peaks  did  not  occur  simultaneously,  even  though  most  current  interpretation 
techniques  for  FWD  data  assume  that  the  maxima  do  occur  simultaneously.  The  farther 
the  measurement  point  was  from  the  loading  point,  the  later  the  peak  deflection  occurred. 
Figure  5-6  illustrates  the  manner  in  which  the  shape  of  the  deflection  basin  changes  with 
time  by  plotting  the  basin  measured  at  6  discrete  times  after  the  impact  occurs.  In  each 
image,  the  thickest  solid  line  represents  the  basing  shape  at  the  stated  time  value.  The 
progressively  thiimer  dashed  lines  represent  the  basin  shapes  at  previous  time  increments. 
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Figure  5-3.  Individual  deflection  responses  due  to  individual  impulse  loads  measured  at 
10  distance  units  from  the  impact  site. 
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Figure  5-4.  Time  history  of  deflection  for  a  point  10  distance  units  from  the  impact  site. 
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Figure  5-5.  a)  Time  history  of  deflection  measured  at  several  radial  distances  from  the  impact  site,  b) 
Deflection  basin  constructed  using  peak  deflections. 

5.1.3  Case  2:  Rolling  Load  Passing  a  Stationary  Point 

This  case  models  a  loaded  wheel  passing  a  stationary  observation  point.  The 
purpose  of  this  C8ise  was  to  develop  the  logic  for  prediction  of  the  time  history  of 
deflection  at  the  observation  point  due  to  the  rolling  load.  Figure  5-7  illustrates  the 
manner  in  which  the  problem  of  a  continuous,  moving  load  was  discretized  as  a  series  of 
impulsive  loads.  Each  point  load,  P,  is  separated  in  time  by  At,  and  in  space  by  Ax.  The 
relationship  between  At  and  Ax  is  a  function  of  the  load  velocity,  V:  Ax  =  VAt.  For  this 
case,  the  load  and  velocity  were  assumed  to  remain  constant  throughout  the  period  of 
interest.  Therefore,  each  impulse  had  equal  magnitude,  I=P*At,  and  was  applied  at  equal 
spacing.  Ax. 

The  same  procedure  used  in  Case  1  to  sum  the  individual  responses  applies  to  this 
case.  However,  in  contrast  to  Case  1,  each  time  history  of  deflection  produced  by  each 
impulsive  load  has  a  different  shape  because  it  is  applied  at  different  distances  from  the 
observation  point.  Figure  5-8  illustrates  the  individual  responses  due  to  each  impulse  for 
two  cases,  one  in  which  the  wheel  passed  directly  over  the  observation  point  (located  at 
x=0,  y=0),  and  the  other,  in  which  the  wheel  traveled  a  parallel  path  offset  by  30  units  in 
the  y-direction  (as  in  Figure  5-7).  The  time  scale  was  adjusted  so  that  at  time  =  0,  the 
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wheel  is  at  x=0,  which  is  directly  over,  or  closest  to,  the  observation  point.  The 
highlighted  individual  deflections  are  caused  by  the  point  load  applied  at  x=0,  when  the 
wheel  is  closest  to  the  observation  point.  As  described  in  Case  1,  the  resulting  overall 
time  history  of  deflection  for  each  case  was  computed  using  superposition,  as  shown  in 
Figure  5-9. 

5.1.4  Case  3:  Deflection  “Snapshot” 

Building  upon  the  results  of  Case  2,  this  case  considered  the  deflections  at  a 
particular  instant  in  time  along  a  line  parallel  to  the  wheel  path.  Such  deflections  could 
be  measured  by  a  series  of  observation  points  aligned  along  the  wheel  path,  if  the 
deflections  were  all  recorded  simultaneously.  In  this  sense,  the  deflections  can  be 
thought  of  as  describing  a  “snapshot”  of  the  deflection  basin  produced  by  a  moving  load. 
The  same  deflection  basin  would  also  be  recorded  by  a  string  of  observation  points  that 
moved  with  the  wheel;  as  such,  this  case  is  more  applicable  to  rolling  wheel 
deflectometers  than  either  of  the  preceding  cases. 

The  snapshot  deflection  basin  was  developed  through  the  following  process. 

First,  a  time  history  deflection  response  was  created  in  the  same  manner  used  in  Case  2. 
Second,  this  time  history  was  shifted  in  time  by  A/  to  represent  the  time  history  of 
deflection  that  would  be  observed  at  (x+Ax,  y).  It  should  be  noted  that  this  is  only  valid 
when  the  material  properties  at  (x+Ax,  y)  are  the  same  as  at  (x,  y).  In  other  words,  an 
observation  point  at  (x+Ax,  y)  would  experience  the  same  time  history  of  deflection  as  an 
observation  point  at  (x,  y),  except  that  the  response  would  start  At  later.  Figure  5-10 
illustrates  this  process  for  several  observation  points  using  the  time  history  for  y=30  from 
Figure  5-9. 
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Figure  5-6.  Deflection  basins  predicted  for  several  discrete  time  intervals. 
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Figure  5-7.  Representation  of  a  rolling  wheel  passing  an  observation  point  as  a  sequence 
of  impulse  point  loads. 
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Figure  5-8.  Individual  impulse  load  deflection  responses  for  case  in  which  wheel  passes  (a)  directly  over 
observation  point  (y=0),  and  (b)  30  distance  units  from  observation  point  (y=30). 
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Time 

Figure  5-9.  Time  history  of  deflection  for  rolling  wheel  passing  directly  over  (y=0)  and 
30  units  (y=30)  from  the  observation  point. 


Figure  5-10.  “Snapshot”  deflection  basin  conceptual  development  illustrating  a  constant  time  plane. 

A  “snapshot”  can  be  created  by  choosing  a  constant  time  value  and  finding  the 
corresponding  deflection  at  each  position.  In  Figure  5-10,  the  dashed-line  vertical  plane 
represents  the  constant  time  slice  for  which  the  snapshot  was  calculated.  Figure  5-11 
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shows  the  resulting  snapshots  for  y=30  and  y=0  when  the  response  calculated  in  Case  2 
(Figure  5-9)  was  shifted  and  cut  as  described  above.  At  the  instant  these  “snapshots” 
were  calculated,  the  wheel  was  directly  over  the  x=0  position.  Note  that  the  peak 
deflection  occurred  slightly  behind  the  position  of  the  wheel  for  both  cases  and  that  the 
peak  deflection  for  y=30  occurred  farther  behind  the  wheel  than  for  the  y=0  case. 


Deflection 


Figure  5-11.  “Snapshot”  deflection  basin  observed  when  wheel  is  directly  over  x  =  0 
position. 

For  this  simplified  case  which  has  constant  velocity  and  load,  the  shape  of  the 
“snapshot”  is  the  reverse  image  of  the  time  history  of  deflection  for  a  stationary 
observation  point,  but  with  an  important  difference  —  the  scale  on  the  x-axis  measures 
distance,  not  time.  When  loads  and  velocities  are  changing  with  time,  this  mirror  image 
effect  will  not  occur  because  the  position  scale  will  not  transform  at  a  constant  rate. 

5.1.5  Case  4:  “Snapshof  ’  3 -Dimensional  Deflection  Basin 

The  “snapshof’  deflection  basin  along  a  single  line  parallel  to  the  wheel  path  was 
developed  in  Case  3.  Case  4  creates  a  three-dimensional  deflection  basin  by  extending 
the  calculations  of  Case  3  to  several  different  parallel  slices  at  different  distances  from  the 
wheel  path.  Plotting  these  snapshots  together  in  three-dimensional  space  illustrates  the 
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shape  of  the  basin  that  travels  with  the  moving  wheel  load  (Figure  5-12).  At  the  instant 
this  basin  was  calculated,  the  load  was  located  at  x=0  and  y=0. 

5.1.6  Case  5:  3-D  Basin  for  Two  Parallel  Rolling  Wheels 

This  final  case  models  two  wheels  traveling  along  parallel  paths  -  i.e.  a  loaded 
axle.  For  this  case,  the  principle  of  superposition  was  applied.  Assuming  the  second 
wheel  applied  the  same  load  as  the  first,  it  would  produce  an  identical  deflection  basin. 
Therefore  a  copy  of  the  single  wheel  basin  can  be  shifted  away  from  the  x-axis  by  a 
distance  equal  to  the  wheel  spacing,  and  superimposed  on  the  original  basin.  Summing 
the  two  basins  creates  the  two-wheel  deflection  basin  illustrated  in  Figure  5-13. 


Figure  5-12.  Deflection  basin  observed  when  a  single  rolling  wheel  is  directly  over  the 
x^O,  y=0  position. 
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Figure  5-13.  Deflection  basin  observed  when  the  axle  connecting  two  rolling 
wheels  is  directly  over  the  x=0  position. 

5.2  Algorithm  Modifications  for  Moving  Deflection  Sensors 

The  previous  section  described  calculating  the  entire  deflection  basin  produced  by 
a  rolling  wheel  at  every  space  increment  along  the  direction  of  travel.  However,  rolling 
weight  deflectometers  have  three  important  differences  from  the  idealized  cases  modeled 
in  Section  5.1. 

First,  rolling  weight  deflectometers  (or  any  other  pavement  testing  device)  only 
measure  deflections  at  several  discrete  points.  Therefore,  only  deflections  at  the 
measurement  sensors  need  to  be  calculated.  This  significantly  reduces  the  required 
computation  time. 

Second,  Case  3  (Section  5.1.4)  assumed  that  the  deflection  sensors  were  stationary 
when  the  “snapshot”  was  measured.  However,  the  deflection  sensors  on  a  rolling  weight 
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deflectometer  move  with  the  load.  This  has  important  implications  for  determining  the 
applied  loads  that  contribute  to  the  deflection  at  the  time  of  measurement. 

Third,  the  previous  section  was  based  on  the  ideal  assumption  that  the  wheel 
rolled  at  constant  velocity  and  applied  a  constant  load.  In  reality,  rolling  weight 
deflectometers  neither  move  at  a  constant  velocity  or  apply  a  constant  load.  Therefore, 
the  model  must  base  its  calculations  on  the  velocity  and  load  recorded  by  the  rolling 
weight  deflectometer  over  the  testing  time. 

The  following  sections  describe  the  algorithm  developed  to  model  a  rolling 
weight  deflectometer.  This  algorithm  was  coded  in  the  computer  program  BASIN. 

5.2.1  BASIN  Algorithm 

A  new  computer  program,  called  BASIN,  was  developed  to  compute  the 
deflections  using  the  load  and  velocity  of  the  rolling  weight  deflectometer.  The  flow 
chart  in  Figure  5-14  illustrates  the  BASIN  algorithm,  which  is  further  described  below. 

To  simplify  the  explanation,  only  a  single  wheel  load  is  considered,  rather  than  a  two¬ 
wheeled  axle. 

First,  data  collected  by  the  rolling  weight  deflectometer  is  loaded  into  the 
program.  This  data  consists  of  the  velocity  and  load  applied  by  the  wheel,  and  the 
deflection  sensor  measurements  at  each  time  increment.  Also,  the  initial  coordinates  of 
the  wheel  (xq,  yo)  and  the  location  of  the  deflection  sensors  relative  to  the  wheel  are 
loaded  into  BASIN.  Because  BASIN  is  a  forward  model,  not  a  backcalculation  tool,  the 
deflection  measurements  are  not  used.  However,  when  a  backcalculation  system  is 
completed,  measured  deflections  could  be  loaded  in  the  algorithm  at  this  point. 

Next,  BASIN  calculates  the  position  of  the  wheel  at  each  time  increment.  When 
the  wheel  location  is  known  at  all  times,  calculation  of  the  position  of  the  defleetion 
sensors  at  all  times  by  adding  or  subtracting  the  relative  offset  of  the  sensors  from  the 
wheel  location  is  straightforward. 
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For  each  deflection  sensor,  there  are  numerous  load  pulses  that  contribute  to  the 
deflection  measured  at  a  given  time.  Each  of  these  pulses  is  applied  at  a  different 
distance  from  the  deflection  sensor;  therefore,  BASIN  next  calculates  all  the  radii  for  the 
applied  load  pulses.  Figure  5-8  illustrated  this  concept. 

Next,  the  time  history  deflection  response  for  each  load  pulse  is  computed  by 
BASIN.  The  response  is  generated  using  the  Green’s  functions  described  in  Chapter  4. 
The  Green’s  function  uses  the  material  properties  of  the  layered  system,  which  are 
supplied  from  an  input  file,  to  calculate  the  deflection  response. 

After  the  deflection  time  history  for  each  load  pulse  is  known,  the  proper 
contribution  to  the  measured  deflection  must  be  selected.  This  is  a  key  step  that  was  not 
necessary  for  the  idealized  cases  modeled  previously.  In  Section  5.1.3  (Case  2),  the 
entire  time  history  of  deflection  due  to  each  applied  load  (at  each  radius)  was  calculated. 
Then  all  time  histories  were  superimposed  to  create  an  overall  time  history  of  deflection 
for  a  single  observation  point.  However,  because  the  time  of  deflection  measurement  is 
known,  only  one  deflection  value  from  each  time  record  is  actually  needed.  Therefore, 
BASIN  selects  the  proper  deflection  contribution  from  each  load  pulse  and  combines  all 
contributions  to  calculate  the  total  deflection. 

Figure  5-15  and  the  following  explanation  help  clarify  this  important  concept. 

When  pulse  Pq  occurs  at  time  Tq,  the  wheel  is  at  position  (%  y^.  The  wheel 
rolls  forward  until  time  TiOi  applying  a  new  load  pulse  at  each  time  step.  At 
time  the  wheel  is  at  location  {xjq,  yj^,  and  the  deflection  sensor  is  at 
(xdefjo,  ydefio).  The  radius  between  Pq  and  the  deflection  sensor  is: 


'lo  =  +  (ydef^o  -  y^^ 

Pq  generates  a  response,  DQ(rjQ),  that  starts  at  time  tQ  and  lasts  through  time  t2o, 
where  the  t  time  scale  is  relative  only  to  each  load  pulse,  and  the  T  scale 
represents  absolute  time.  Only  the  deflection  at  time  tjg  contributes  to  the 
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total  deflection  measured  by  the  deflection  sensor.  Similarly,  P;  creates  a 
response  and  contributes  Di(rg,tg)  at  time  Tjo-  Obviously,  only  pulses 
occurring  prior  to  the  measurement  time  can  influence  the  measured 
deflection,  so  the  last  pulse  that  must  be  considered  is  Pp.  Therefore  the  total 
deflection  measured  by  the  sensor  at  time  Tio  is: 


(^0  ^0  (^10  5  ho  A  ^9 » h  A  (^8  5  4  )*■•  •  •  •'*' A  ^2  ’  h  A  (h  ’  h 


Figure  5-15.  Method  for  calculating  deflections  at  sensor  location  Dj  at  time  Tjo-  For  simplification  of  the 
illustration,  haversine  load  pulses  are  represented  by  vertical  arrows. 
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After  the  last  load  pulse  is  considered,  BASIN  repeats  the  identical  process  for 
the  next  deflection  sensor.  An  entirely  new  set  of  radii  and  responses  must  be  calculated 
because  each  sensor  location  is  different.  When  the  total  deflection  for  the  final  sensor 
has  been  calculated,  BASIN  moves  to  the  next  set  of  measured  deflections  and  calculates 
the  basin  measured  at  the  next  time  increment.  For  this  thesis,  the  load,  velocity,  and 
material  properties  were  assumed  to  be  constant  throughout.  Therefore,  only  one  basin 
was  calculated  because  all  others  would  be  identical. 

5.2.2  Avetime  and  Timehack 

To  limit  the  number  of  required  computations,  two  questions  must  be  answered 
regarding  the  load  pulses.  First,  how  far  apart  (in  time  and  space)  should  the  load  pulses 
occur?  Second,  how  far  back  (in  time  and  space)  should  loads  be  considered  to  have  an 
influence  on  the  measured  deflection?  Ideally  loads  would  be  placed  infinitesimally 
close  and  all  previous  loads  would  be  considered.  However,  this  is  computationally 
impractical.  Two  parameters  were  created  to  quantify  these  two  issues:  avetime  and 
timeback. 

Avetime  represents  the  time  increment  between  each  load  pulse.  For  the  idealized 
cases  analyzed  in  Section  5.1,  point  loads  were  assumed  to  be  applied  at  every  time  (and 
space)  increment.  When  the  transition  to  the  real  Green’s  function  solution  was  made, 
haversine  load  pulses  were  applied  such  that  they  overlapped  to  produce  a  constant  line 
load.  Section  4.2.4.4  first  addressed  this  issue  and  showed  that  different  combinations  of 
spacing  and  pulse  amplitude  could  produce  the  same  line  load  in  time.  The  best  spacing 
between  load  pulses  was  only  suitably  determined  through  a  parametric  study.  Section 
7.2.5  describes  the  manner  in  which  the  predicted  deflection  basin  changes  as  the  spacing 
is  varied. 

Timeback  represents  the  period  of  time  over  which  load  pulses  are  considered  to 
influence  the  calculated  deflections.  This  has  a  significant  influence  on  the  computation 
time  because  the  number  of  load  pulses  increase  in  direct  proportion  to  timeback.  By 
examining  time  history  deflections  predicted  by  FWDDROP,  it  becomes  clear  that  for  an 


70 


average  load  pulse  duration  of  30  msec,  all  responses  are  negligible  after  200  msec. 
Section  7.2.4  examines  this  parameter  in  greater  detail  to  determine  its  influence  on  the 
deflection  basin  shape. 

5.2.3  Deflection  Basin  Assembly 

When  the  deflection  at  each  measurement  location  is  known,  the  data  can  be 
assembled  and  plotted  to  show  the  approximate  shape  of  the  deflection  basin.  Two-  or 
three-dimensional  images  can  be  created,  depending  on  the  location  of  the  deflection 
sensors.  Figure  5-16  illustrates  a  longitudinal  deflection  basin  computed  by  BASIN.  The 
data  point  markers  show  the  location  of  the  deflection  sensors  on  a  line  directly  under  the 
wheel  path.  Note  that  more  sensors  were  located  at  points  of  high  curvature  to  more 
accurately  define  the  basin. 

Using  the  method  described  for  the  Case  4  idealized  model  (Section  5.1.5),  three- 
dimensional  basins  can  be  constructed  using  BASIN.  Figure  5-17  illustrates  such  a  basin. 


Distance  [m] 


Figure  5-16.  Longitudinal  deflection  basin  with  deflection  sensors  located  at  each  data  point  marker. 
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Figure  5-17.  Single  wheel  deflection  basin  computed  using  BASIN. 

5.2.4  3-D  Basin  for  Two  Parallel  Rolling  Wheels 

The  deflection  basin  created  by  two  parallel  rolling  wheels,  each  with  identical 
load,  can  be  constructed  using  superposition.  As  described  in  Case  5  (Section  5.1.6),  a 
two-wheel  basin  is  created  when  a  copy  of  the  one- wheel  basin  is  shifted  in  the  transverse 
direction  by  a  distance  equal  to  the  wheel  spacing,  and  then  superimposed  on  the  original 
basin.  Figure  5-18  illustrates  the  result  of  this  process. 
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Figure  5-18.  Double  wheel  deflection  basin  computed  using  BASIN. 


5.4  BASIN  Code  Details 

The  idealized  cases  examined  in  Section  5.1  were  programmed  using  MATLAB 
because  of  the  user-friendly  syntax  and  enhanced  graphics  capabilities.  However,  when 
the  complex  Green’s  function  was  used,  the  program  was  converted  to  MS-FORTRAN. 
The  entire  code  is  listed  in  Appendix  A. 


Chapter  6 

PAVEMENT  DEFLECTION  MODEL  VERIFICATION 


To  test  the  accuracy  of  BASIN,  the  program  was  compared  to  other  pavement 
models.  Data  from  dynamic  moving  wheel  models  was  not  available  for  comparison, 
therefore  BASIN  was  checked  against  static  pavement  deflection  models.  To  make  this 
comparison,  BASIN  was  modified  to  calculate  the  deflections  that  would  be  produced  by 
a  FWD  test.  Using  moduli  that  were  backcalculated  from  static  models,  the  deflections 
from  the  new  program,  FWDDROP,  were  compared  to  deflections  measured  by  FWD 
tests. 


6.1  FWDDROP  Program 

To  create  a  program  that  would  predict  the  pavement  deflections  produced  by  a 
falling  weight  deflectometer,  the  BASIN  program  was  modified  to  compute  the  time 
history  of  deflection  at  a  several  fixed  locations  resulting  from  a  transient  load  applied 
over  a  known  area  at  a  known  location.  Several  distinct,  yet  simple,  changes  were 
required  to  accomplish  the  modification.  These  changes  resulted  in  a  new  computer 
program,  FWDDROP.  The  steps  required  to  develop  FWDDROP  are  described  in  the 
following  paragraphs. 

First,  basin’s  data  input  module  was  scaled  down  because  the  input  data  required 
for  an  FWD  model  is  less  comprehensive  than  that  required  for  the  RWD  model.  For 
FWDDROP,  only  the  following  items  need  to  be  known.  (Units  for  the  program  can  be 
either  metric  or  US,  as  long  as  they  are  consistent.) 

Pavement  system  parameters  -  variables  needed  to  define  the  geometry  and 
material  properties  of  each  layer  of  the  pavement  profile  (thickness,  shear  wave 
velocity,  mass  density,  Poisson’s  ratio,  and  damping  ratio). 

Load  -  the  peak  load  applied  by  the  FWD. 

Load  Duration  -  the  period  of  time  over  which  the  impulse  load  is  applied. 
Load  Area  -  the  size  of  the  FWD  loading  plate. 

Sensor  Locations  -  the  locations  of  the  FWD  deflection  sensors  (geophones). 
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Second,  the  BASIN  modules  that  determine  the  location  of  the  wheel  and  deflection 
sensors  were  removed  because  the  load  plate  and  the  sensors  are  stationary  throughout  the 
test.  Third,  the  program  was  changed  so  that  only  one  load  pulse  is  applied  to  the 
pavement,  rather  than  the  multiple  pulses  used  to  model  a  moving  load.  This  modification 
required  removal  of  the  loop  that  calculated  the  response  from  each  pulse.  The  coding  in 
BASIN  that  characterizes  the  load  as  a  haversine  was  retained  in  FWDDROP;  however, 
any  load  shape  (or  an  actual  load  time  history)  could  be  used  in  FWDDROP.  Finally,  the 
output  data  module  was  modified  to  save  the  entire  time  history  of  deflection  at  each 
sensor  location. 

With  the  modifications  complete,  FWDDROP  was  a  fully  functional  dynamic 
linear-elastic  pavement  deflection  model.  Figure  6-1  illustrates  the  deflection  time 
histories  produced  by  FWDDROP  and  shows  how  a  “static”  deflection  basin  can  be 
produced  by  selecting  the  peak  value  from  each  deflection  time  history,  then  plotting 
deflection  as  a  function  of  distance  from  the  center  of  the  load  plate. 


Figure  6-1.  Typical  FWD  deflection  time  histories  and  construction  of  “static”  deflection  basin. 
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6.2  Comparison  to  SHRP  Test  Sections 

Deflections  calculated  by  FWDDROP  were  compared  to  deflections  reported  by 
Meier  (1995)  for  the  two  Strategic  Highway  Research  Program  (SHRP)  pavement 
sections  illustrated  in  Figure  6-2.  Meier  used  three  static,  layered  elastic  pavement 
backcalculation  programs  to  determine  the  moduli  of  the  SHRP  pavement  layers,  as  listed 
in  Table  6-1 .  The  first  set  of  moduli  was  backcalculated  using  Meier’s  artificial  neural 
network  algorithm  which  he  coded  in  an  extremely  fast  basin-matching  program.  The 
second  set  of  moduli  were  determined  using  Modulus  4.0,  SHRP’s  database-searching, 
basin-matching  program  (Rada,  Richter,  and  Jordahl,  1994).  The  last  set  of  moduli  were 
backcalculated  by  WESDEF,  an  iterative  program  that  uses  a  layered  elastic  pavement 
model  that  allows  variable  frictional  properties  on  layer  interfaces  (Van  Cauwelaert  et  al., 
1989).  Considerable  scatter  (5%)  exists  in  the  AC  layer  moduli  for  both  sections,  even 
though  all  analyses  were  performed  with  static  load  models. 


SHRP  Section  A 

Asphalt  Concrete 

4.95  in  =0.126  m 

Crushed  Limestone  Base 

13.4  in  =0.340 

m 

Soil/Aggregate  Subbase 

12.0  in  =0.305  m 

Sand  Subgrade 

00 

SHRP  Section  B 

Asphalt  Concrete 

4.20  in  =0.1 07  m 

Uncrushed  Gravel  Base 

5.00  in  =0. 127  m 

Sand  Subgrade 

00 

Figure  6-2.  Pavement  structures  for  SHRP  test  sections  A  and  B  (after  Meier,  1995). 

Because  all  the  parameters  of  the  pavement  system  and  the  loading  time  history 
were  not  available,  some  values  were  estimated  from  the  best  available  information. 
Estimated  values  of  mass  density,  Poisson’s  ratio,  and  damping  ratio  for  each  layer  of  the 
pavement  system  are  listed  in  Table  6-2.  Mass  densities  were  based  on  suggested  values 
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reported  by  Yoder  and  Witczak  (1975).  Poisson’s  ratios  were  taken  as  those  assumed 
by  Meier  (1995)  for  his  backcalculations.  The  damping  ratio  was  assumed  to  be  5%  for 
all  materials.  Sensitivity  calculations  showed  that  errors  in  the  assumed  density  or 
damping  ratio  had  almost  no  effect  on  the  calculated  deflection  basins. 


Table  6-1.  Moduli  backcalculated  from  measured  deflection  basins 
(Meier,  1995).  [1  ksi  =  6.895  MPa] 


Backcalculated  Moduli  (ksi) 

SHRP 

Pavement 

Neural 

Modulus 

Section 

Layer 

Network 

4.0 

WESDEF 

A 

Asphalt 

1294 

1250 

1317 

Base* 

42 

41 

42 

Subgrade 

32 

30 

31 

B 

Asphalt 

855 

921 

918 

Base 

53 

56 

46 

Subgrade 

27 

27 

27 

‘Combination  of  crushed  limestone  base  and  soil/aggregate  subbase 


Table  6-2.  Assumed  layer  properties  for  SHRP  test  sections. 
[1  Ib/ft^  =  16.018  kg/m^] _ 


SHRP 

Section 

Pavement 

Layer 

Density 

[Ib/fP] 

Poisson's 

Ratio 

Damping 

Ratio 

A&B 

Asphalt 

150.00 

0.325 

0.05 

Base 

130.00 

0.35 

0.05 

Subgrade 

120.00 

0.35 

0.05 

The  load  pulse  was  modeled  as  a  haversine  with  duration  of  26.5  msec  as  suggested 
by  Meier  (1995).  Meier  found  that  differences  between  the  shape  of  an  actual  load  pulse 
and  a  haversine  required  that  the  peak  of  the  haversine  be  increased  by  a  factor  of  1 .09  to 
achieve  the  same  total  impulse.  Therefore,  the  peak  of  the  haversine  applied  in 
FWDDROP  was  9%  greater  than  the  actual  metisured  load. 
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Deflection  basins  were  calculated  using  each  set  of  moduli  for  both  pavement 
sections,  and  the  error  between  the  measured  and  calculated  deflections  was  determined. 
Table  6-3  presents  the  results  of  the  calculations.  Figure  6-3  shows  the  measured 
deflection  basins  and  the  basins  predicted  by  FWDDROP  using  Meier’s  neural  network 
backcalculated  moduli.  The  agreement  between  the  deflection  basins  is  very  good, 
particularly  when  all  of  the  assumptions  that  went  into  the  dynamic  model  are  considered. 
It  should  be  noted  that  the  moduli  used  by  FWDDROP  were  backcalculated  using  static 
pavement  deflection  models,  and  the  true  dynamic  modulus  of  each  layer  is  unknown. 


Table  6-3.  Comparison  of  measured  and  calculated  FWD  deflections  for  SHRP  test  sections. 
[1  mil  =  25.4  pm] _ 


Section 

Layer  Moduli 

Model 

Measured 

1 

1 

1 

I 

1 

Load  [lbs] 

0 

60 

A 

9.47 

7.43 

6.31 

3.82 

2.48 

1.24 

10006 

Neural  Network 

Calculated 

9.28 

7.47 

6.26 

4.78 

3.72 

2.47 

1.48 

error  (%) 

-2.0 

0.6 

-0.8 

-2.4 

-2.7 

-0.3 

19.0 

Modulus  4.0 

Calculated 

9.63 

7.76 

6.51 

4.98 

3.89 

2.61 

1.56 

error  (%) 

1.6 

4.4 

3.1 

1.7 

1.9 

5.1 

26.0 

WESDEF 

Calculated 

9.31 

7.52 

6.32 

4.85 

3.78 

2.53 

1.52 

error  (%) 

-1.7 

1.2 

0.1 

-1.1 

-0.9 

2.1 

22.2 

WMWIHilliMliM 

HHHHB 

B 

Measured 

12.21 

9.43 

7.65 

5.72 

4.32 

2.62 

1.36 

9596 

Neural  Network 

Calculated 

12.46 

9.64 

7.89 

5.92 

4.57 

2.99 

1.72 

error  (%) 

2.0 

2.2 

3.2 

3.4 

5.7 

13.9 

26.7 

Modulus  4.0 

Calculated 

12.14 

9.47 

7.80 

5.89 

4.56 

2.99 

1.72 

error  (%) 

-0.6 

0.4 

1.9 

3.0 

5.5 

14.0 

26.6 

WESDEF 

Calculated 

12.56 

9.77 

7.99 

4.59 

2.99 

1.72 

error  (%) 

2.9 

3.6 

4.5 

■ 

6.2 

13.9 

26.8 

6.3  Comparison  to  PACCAR  Test  Section 

An  additional  pavement  section  was  used  to  further  evaluate  the  accuracy  of 
FWDDROP.  A  heavily  instrumented  asphalt  concrete  pavement  section  located  at 
PACCAR  Inc.’s  Technical  Center  in  Moimt  Vernon,  WA,  was  extensively  tested  and 
analyzed  by  Winters  (1993).  Data  collected  from  130  FWD  (WSDOT’s  Dynatest  8000) 
tests  at  61  locations  on  a  12.2  m  (40  ft)  pavement  section  were  analyzed  using 
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EVERCALC  3.3  (a  static,  layered  elastic  backcalculation  program)  to  determine  layer 
moduli. 


SHRP  Section  A 


Figure  6-3.  Measured  and  calculated  deflection  basins  for  SHRP  sections  A  and  B. 
[1  mil  =  25.4  nm] 


Measured  deflections  from  three  FWD  tests  were  compared  with  FWDDROP 
calculations.  All  three  tests  were  conducted  at  the  site  identified  as  Core  4,  where  layer 
thicknesses  had  been  determined  from  a  drilled  core,  illustrated  in  Figure  6-4.  Winters 
backcalculated  the  moduli  of  the  AC,  base,  and  subgrade  for  each  test,  as  shown  in  Table 
6-4.  One  significant  difference  between  the  PACCAR  profile  and  the  SHRP  sites  is  that 
the  material  below  the  groundwater  table  was  modeled  as  a  stiff  layer.  Using 
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EVERCALC,  Winters  determined  that  the  modulus  of  the  stiff  layer  was  40  ksi  (275.8 
MPa).  Other  researchers  agree  that  a  fully  saturated  soil  layer  should  be  modeled  as  a 
stiff  layer  for  backcalculation  purposes  (Newcomb  et  ah,  1995). 


PACCAR  Core  4 


Asphalt  Concrete 

5.4”=0.137m 

Crushed  Stone  Base 

13.0”=0.33m 

Sandy  Clay  Subgrade 

46.1”=1.171  m 

Saturated  Stiff  Layer 

00 

Figure  6-4.  Pavement  structures  for 
PACCAR  test  site  Core  4. 


Table  6-4.  Moduli  backcalculated  from  measured  deflection  basins. 
[1  ksi  =  6.895  MPa] 


PACCAR 

FWD 

Backcalculated  Moduli  (ksi) 

Cored 

Load  (lb) 

Asphalt 

Base 

Subgrade 

1 

14039 

654.794 

11.42 

9.542 

2 

10631 

678.106 

11.634 

10.14 

3 

5045 

731.217 

10.688 

13.185 

Winters  reported  the  assumed  Poisson’s  ratios  used  for  his  backcalculations,  but  the 
densities  and  damping  ratios  were  not  known.  These  values  were  again  assumed  as 
described  in  Section  6.2,  and  are  listed  in  Table  6-5. 

The  load  pulse  was  modeled  as  a  30  msec  haversine,  which  is  consistent  with  the 
duration  reported  by  Chatti  et  al.  (1995)  for  WSDOT’s  Dynatest  8000.  The  reported  load 


pulse  for  WSDOT’s  FWD  very  closely  resembled  a  haversine,  therefore  no 
adjustments  to  match  the  total  impulse  were  required. 
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Table  6-5.  Assumed  layer  properties  for  PACCAR  Core  4. 
[1  Ib/ft^  =  16.018  kgV] _ 


PACCAR 

Core  4 

Pavement 

Layer 

Density 

[ib/fP] 

Poisson's 

Ratio 

Damping 

Ratio 

1.2,3 

Asphalt 

150.00 

0.350 

0.05 

Base 

130.00 

0.40 

0.05 

Subgrade 

110.00 

0.45 

0.05 

Stiff 

120.00 

0.35 

0.05 

Table  6-6  lists  the  measured  and  predicted  deflections  and  errors,  and  Figure  6-5 
illustrates  the  deflection  basins.  The  predicted  deflections  do  not  agree  with  the  measured 
deflections  as  closely  as  they  did  for  the  SHRP  pavement  sections;  however,  the  shape  of 
the  calculated  basin  is  virtually  parallel  to  that  of  the  measured  basin.  It  would  appear 
that  the  basins  would  match  very  closely  if  the  peak  load  was  reduced.  However,  Figure 
6-6  illustrates  that  an  8%  reduction  in  load  changes  the  deflections  differently  at  each 
distance  from  the  load;  therefore,  the  calculated  basin  loses  it’s  parallel  match  to  the 
measured  basin.  The  discrepancies  between  the  basins  can  be  possibly  attributed  to  the 
dynamic  interaction  of  the  stiff  layer.  As  described  in  Chapter  3,  stiff  layers  at  shallow 
depths  can  profoundly  affect  the  measured  deflections,  which  can  result  in  significant 
errors  from  backcalculation  programs  that  use  static  loading  models. 

6.4  Modification  of  Deflection  Time  Histories 

During  analysis  of  the  deflection  time  histories,  a  problem  was  encountered  with  the 
deflections  at  the  very  beginning  and  end  of  the  time  record.  Figure  6-7a  indicates  that 
the  initial  pavement  deflection  is  positive  at  the  time  of  FWD  impact.  This  implies  that 
the  pavement  has  moved  before  loading  has  been  applied.  Also,  the  deflections  at  the 
very  end  of  the  time  record  become  positive,  again  implying  that  the  pavement  begins  to 
deflect  again,  long  after  the  deflections  due  to  the  loading  have  ceased.  Presumably, 


these  spurious  deflections  are  due  to  a  basic  assumption  in  a  Fourier  transformation  — 
the  motion  transformed  into  the  fi'equency  domain  is  periodic  (i.e.,  the  motion  begins 
again  immediately  after  the  last  point  of  the  time  record).  It  is  also  possible  that  the 
spurious  deflections  may  result  in  part  from  the  Green’s  function  assumption  that  the 
response  varies  linearly  between  the  top  and  the  bottom  of  each  sublayer. 


Table  6-6.  Comparison  of  measured  and  calculated  FWD  deflections  for  PACCAR  Core  4. 
[1  mil  =  25.4  pm] _ 


Deflections  [mils]  at  radius  [in] 

Test 

Load  [lbs] 

0 

8 

12 

24 

36 

48 

1 

14039 

Measured 

29.20 

24.06 

20.58 

11.66 

6.04 

3.17 

Calculated 

31.18 

26.42 

22.97 

14.08 

8.37 

5.16 

error  (%) 

6.8 

9.8 

11.6 

20.8 

38.6 

62.8 

• 

2 

10631 

Measured 

21.39 

17.57 

15.04 

8.52 

4.41 

2.33 

Calculated 

22.64 

19.21 

16.67 

10.18 

6.03 

3.71 

error  (%) 

5.9 

9.4 

10.9 

19.5 

36.8 

59.4 

3 

5045 

Measured 

9.29 

7.62 
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Input  parameters  were  systematically  varied  to  determine  which  affected  this 
problem.  The  thickness  of  the  layers  of  the  pavement  system  appeared  to  have  the 
greatest  influence,  but  the  spurious  deflections  could  not  be  completely  removed  by 
reducing  the  thickness  of  individual  layers. 


Deflection  [mils] 


Figure  6-5.  Measured  and  calculated  deflection  basins  for  PACCAR  Core  4 
FWD  tests.  [1  mil  =  25.4  fim] 


Deflection  [mils]  Deflection  [mils]  Deflection  [mils] 


Figure  6-6.  Measured  and  calculated  deflection  basins  for  PACCAR  Core  4 
FWD  tests  when  peak  load  is  reduced  by  8%.  [1  mil  =  25.4  pm] 
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A  correction  procedure  was  developed  to  eliminate  the  computed  spurious 
deflections.  A  simple  “tripwire”  code  loop  was  added  to  FWDDROP  (and  BASIN)  to 
process  the  calculated  time  history  of  deflection.  From  the  beginning  of  the  time  record, 
the  deflection  is  set  to  0.0  until  the  first  negative  deflection  occurs.  From  that  time  until 
the  middle  of  the  time  record,  the  deflections  are  not  modified.  By  the  middle  of  the  time 
record,  pavement  deflections  are  negligible,  so  the  deflections  were  set  to  0.0  for  the 
remainder  of  the  loading  period.  Figure  6-7b  illustrates  the  modified  time  history  of 
deflection  calculated  for  the  same  sensor  used  in  Figure  6-7a.  The  peak  deflections 
calculated  by  FWDDROP  remain  absolutely  unchanged  by  this  modification.  Also, 
comparisons  between  deflection  basins  calculated  by  BASIN  using  the  modified  and 
unmodified  deflection  records  revealed  only  very  minor  differences. 

Deflection  [mils]  Deflection  [mils] 


Figure  6-7.  a)  Illustration  of  negative  deflections  at  the  beginning  and  end  of  the  deflection  time  history 
calculated  by  FWDDROP.  b)  Modified  deflection  time  history. 


Chapter  7 

PARAMETRIC  STUDY 


The  sensitivity  of  the  calculated  deflection  basin  to  a  variety  of  input  parameters 
was  investigated  by  a  detailed  parametric  study.  The  parametric  study  involved 
systematic  perturbation  of  each  input  variable  and  evaluation  of  the  resulting  changes  in 
calculated  deflections.  This  chapter  describes  the  plan,  execution,  and  results  of  the 
study. 

7.1  Parametric  Study  Plan 

A  detailed,  logical  plan  of  action  for  the  parametric  study  was  prepared  before 
calculations  began.  First,  a  standard  pavement  section  and  standard  loading  were  selected 
for  use  throughout  the  study.  Next,  the  BASIN  inputs  were  grouped  into  three  categories 
of  variables  according  to  their  role  in  the  calculations;  algorithm  variables,  pavement 
system  variables,  and  loading  variables.  Last,  a  method  was  developed  for  quantitatively 
describing  and  comparing  the  calculated  deflection  basins.  Each  of  these  aspects  of  the 
plan  are  described  in  the  following  sections. 

7.1.1  Standard  Pavement  Section  and  Standard  Loading 

A  pavement  section  was  selected  to  establish  a  standard  deflection  basin  that  would 
serve  as  a  common  starting  point  for  all  sensitivity  analyses.  To  keep  the  basin 
consistent,  a  standard  loading  scheme  was  also  designated. 

A  12.2  m  (40  ft)  section  of  pavement  at  PACCAR  Inc.’s  Technical  Center  in  Mt. 
Vernon,  WA,  was  selected  as  the  standard  pavement  section  for  two  primary  reasons. 
First,  it  is  a  very  typical  heavy-traffic  flexible  pavement  structure  —  the  thickness  and 
properties  of  the  asphalt,  base,  and  subgrade  layers  are  very  typical  of  those  commonly 
used  throughout  Washington  state.  Second,  the  sections  had  been  extensively  analyzed 
by  Winters  (1993);  consequently,  the  layer  thicknesses  and  moduli  were  well-established 
from  core  samples  and  backcalculations  of  over  one  hundred  FWD  tests. 
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Figure  7-1  illustrates  the  profile  of  the  PACCAR  test  section.  The  pavement 
structure  consists  of  137  mm  (5.4  in)  of  asphalt  on  330  mm  (13.0  in)  of  crushed  stone 
base.  The  subgrade  material  is  a  medium  stiff  sandy  clay.  The  groundwater  table  was 
measured  at  the  time  of  FWD  testing  and  found  at  an  average  depth  of  1.637  m  (64.4  in) 
below  the  pavement  surface.  The  modulus  of  each  layer  was  determined  from  the 
average  of  130  FWD  analyses  using  EVERCALC  3.3.  As  described  in  Chapter  6, 
Winters  modeled  the  saturated  soil  below  the  ground  water  table  as  a  stiff  layer  with  a 
modulus  of  275.8  MPa  (40  ksi).  The  density  of  each  layer  was  assumed  based  on 
suggested  values,  and  the  Poisson’s  ratios  were  set  equal  to  the  values  used  by  Winters 
(1993)  for  backcalculation.  A  damping  ratio  of  5%  was  assigned  to  each  layer. 

In  addition  to  the  standard  pavement  section,  a  standard  load  was  also  selected.  The 
standard  wheel  load  was  set  equal  to  the  load  applied  by  the  left  front  wheel  of  a 
Peterbuilt  truck,  model  PB359.  The  standard  tire  pressure  was  selected  at  90  psi,  which 
is  a  typical  operating  pressure  for  large  trucks  in  the  United  States.  This  combination  of 
load  and  tire  pressure  results  in  a  standard  loaded  radius  of  0.1 12  m  (4.4  in).  As 
explained  in  Chapter  4,  the  standard  load  pulse  has  a  haversine  shape  with  a  duration  of 
30  msec.  This  duration  corresponds  to  a  standard  travel  speed  of  50  km/hr  =  13.89  m/sec 
(31  mph). 

7.1.2  Structure 

The  input  variables  of  BASIN  were  grouped  into  three  categories:  algorithm 
variables,  pavement  system  variables,  and  loading  variables. 

Algorithm  variables  result  from  the  concepts  and  mathematics  used  in  BASIN. 

They  include  the  time  increment,  frequency  interval,  frequency  range,  layer 
discretization,  spacing  between  load  pulses,  and  the  number  of  pulses  behind  the  wheel 
that  are  eonsidered  in  the  deflection  calculations.  These  parameters  were  evaluated  first 
beeause  they  influence  the  speed  and  accuracy  of  the  deflection  calculations.  The 
algorithm  variables  were  adjusted  to  find  values  that  produce  high  accuracy  with  minimal 
computational  effort. 


Figure  7-1.  Standard  Pavement  Section. 
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Pavement  system  variables  describe  the  pavement  structure  that  supports  the 
wheel  load.  They  include  the  thickness,  modulus,  density,  Poisson’s  ratio,  and  damping 
ratio  of  each  layer.  In  typical  static  pavement  analyses,  all  except  the  layer  moduli  are 
assumed.  However,  because  BASIN’ s  dynamic  analysis  may  be  more  sensitive  to  the 
other  parameters,  each  of  these  variables  were  studied  separately. 

Loading  variables  characterize  the  magnitude  and  type  of  load  applied  to  the 
pavement.  They  include  the  magnitude,  pulse  dvuation,  vehicle  speed,  and  tire  pressure. 
These  parameters  are  purely  functions  of  the  loading  mechanism;  however,  it  is  important 
to  understand  their  effects  on  the  basin  shape  because  the  sensitivity  of  the  deflection 
basin  to  the  loading  variables  may  help  determine  the  optimum  loading  conditions  and 
the  best  locations  for  deflection  sensors. 

7.1.3  Data  Presentation  and  Statistical  Parameters 

Because  the  sensitivity  of  deflection  basins  to  parameter  variations  is  not  easily 
described  in  qualitative  terms,  a  quantitative  approach  was  developed.  For  this  purpose, 
the  shape  of  the  deflection  basin  was  described  by  a  series  of  coefficients  related  to  the 
moments  of  the  deflection  basin.  Description  of  distributed  parameters  by  moments  is 
frequently  used,  most  commonly  in  the  area  of  statistics.  If  the  deflection  basin  is 
considered  analogous  to  a  probability  density  function,  the  mean,  standard  deviation, 
coefficient  of  skewness,  and  coefficient  of  kurtosis  can  be  calculated.  These  parameters 
are  most  easily  described  and  understood  by  referring  to  Figure  7-2. 


Figure  7-2.  Example  function,  z(x),  represented  by  N  discrete  points. 
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The  area  under  the  function  can  be  calculated  using  the  following  formula: 

Area  =  ^z,.Ay,. 

/=i 

Next,  a  function, analogous  to  the  probability  density  function,  can  be  calculated  for 
all  Xj  to  show  the  likelihood  that  a  randomly  sampled  value  of  z  will  be  at  x. 

.  _ 

Area 

After  fi  is  determined  for  all  N  terms,  several  statistical  parameters  can  be  calculated. 


Mean,  x  :  the  mean,  average,  or  expected  value  of  the  function  is  calculated  as: 

x  =  ixJi 
1=1 

Figure  7-3  shows  how  the  mean  for  two  identical  functions  can  be  different, 
depending  on  the  location  of  the  distribution  on  the  x  axis. 


Figure  7-3.  Two  functions  with  identical  shapes  but 
different  mean  values. 


Standard  Deviation,  a:  the  standard  deviation  is  a  measure  of  the  spread  of  the 
function  about  the  mean  and  is  related  to  the  variance,  V,  or  second  moment  about 
the  mean,  The  standard  deviation  is  equal  to  the  square  root  of  the  variance, 
and  is  calculated  in  the  following  manner. 
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/=1 

a  =VV 

Because  the  standard  deviation  h2is  the  same  units  as  x,  its  actual  value  is  only 
meaningful  when  the  units  of  x  are  known.  Figure  7-4  illustrates  two  functions 
centered  about  the  same  mean,  but  with  different  distributions  —  one  broad  (large 
ct),  one  narrow  (small  a). 


Figure  7-4.  Two  functions  with  identical  means  but 
different  standard  deviations. 


Coefficient  of  Skewness,  C^:  The  coefficient  of  skewness  provides  a  quantitative 
measure  of  the  symmetry  of  the  function  about  its  mean  value.  The  coefficient  of 
skewness  is  calculated  from  the  third  moment  about  the  mean,  ^,5.  To  render  the 
coefficient  of  skewness  independent  of  scale,  the  third  moment  is  divided  by  cube 
of  the  standard  deviation. 

N 

t=i 


Skewness  can  be  either  positive,  negative,  or  zero.  If  more  than  half  of  the 
function  is  greater  than  the  mean,  the  function  is  negatively  skewed.  If  the 
majority  of  the  function  is  less  than  the  mean,  the  function  is  positively  skewed 
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(Walker  and  Lev,  1969),  as  illustrated  in  Figure  7-5.  A  symmetric  function  has 
zero  skewness. 


Figure  7-5.  Two  functions  with  opposite  skewness,  one 
positive  and  the  other  negative. 


Coefficient  of  Kurtosis,  C4:  Kurtosis  is  a  measure  of  the  peakedness  of  a  function 
calculated  from  the  fourth  moment  about  the  mean,  The  dimensionless 
coefficient  of  kurtosis  is  equal  to  the  fourth  moment  divided  by  the  fourth  power 
of  the  standard  deviation. 

/=i 


Kurtosis  describes  the  arch  of  the  distribution,  and  is  usually  compared  to  a 
normal  distribution  with  a  variance  equal  to  the  variance  of  the  distribution  of 
interest.  The  coefficient  of  kurtosis  of  a  normal  distribution  is  always  equal  to  3. 
Therefore,  if  the  peak  of  the  studied  function  is  higher  than  the  normal 
distribution  peak,  C4  is  greater  than  3 .  If  the  function  is  flatter  on  top  and  with  a 
lower  peak  than  a  normal  distribution,  C4  is  less  than  3  (Sachs,  1982).  Kurtosis  is 
illustrated  in  Figure  7-6. 
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.C4>3 


Figure  7-6.  A  normal  distribution  (C4  =  3)  and  two  functions  with 
greater  and  lesser  kurtosis. 

The  results  of  each  parameter  study  are  illustrated  with  the  sequence  of  figures 
shown  in  Figure  7-7.  Figure  7-7a  shows  the  longitudinal  deflection  basin  for  each  of  the 
parameters  of  interest.  Longitudinal,  rather  than  transverse,  basins  were  computed 
throughout  this  study  because  they  were  expected  to  reveal  more  information  about 
dynamic  deflection  behavior.  In  all  deflection  basin  figures,  the  loaded  wheel  is  located 
at  the  zero  position  coordinate  and  moving  from  left  to  right.  Figure  7-7b  shows  how  the 
peak  deflection  varies  for  each  parameter  value.  Figure  7-7c  illustrates  the  mean  and 
standard  deviation  of  each  deflection  basin.  The  coefficient  of  skewness  is  plotted  in 
Figure  7-7d,  and  the  coefficient  of  kurtosis  is  illustrated  in  Figure  7-7e.  In  Figure  7-7a, 
the  standard  deflection  basin  is  always  printed  using  a  bold  line.  In  Figures  7-7b  through 
7-7e,  the  standard  value  is  designated  by  a  solid  data  point  marker. 

7.2  Algorithm  Variables 

The  algorithm  variables  are  the  basic  limits  used  when  BASIN  computes  the 
deflection  basins.  They  can  have  a  profound  impact  on  the  accuracy  of  the  results  and  the 
computation  time  required  to  achieve  them.  These  parameters  were  investigated  first  to 
find  the  values  that  maximized  accuracy  while  minimizing  computation  time. 


Longitudinal  Deflection  Basins 
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Peak  Deflection 


Mean  and  Std.  Deviation 
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Figure  7-7.  Typical  data  presentation  format. 

7.2.1  Layer  Discretization 

The  Green’s  function  calculation  scheme  developed  by  Kausel  and  applied  in 
PUNCH  relies  on  the  individual  layers  being  sufficiently  thin  that  a  linear  response 
within  a  layer  closely  approximates  the  actual  response.  Thus,  if  layers  are  too  thick,  the 
response  may  be  in  error. 

To  determine  the  best  layering  scheme  for  the  standard  pavement  section,  the 
number  of  sublayers  in  each  of  the  four  layers  was  initially  set  to  one.  Then  the  number 
of  sublayers  was  systematically  increased  until  there  were  a  total  of  40  sublayers  in  the 
entire  pavement  system.  A  total  of  six  layering  schemes  were  tested. 

As  the  number  of  sublayers  was  increased,  the  basin  converged  to  a  constant  shape 
(Figure  7-8a).  The  basin  shape  appeared  to  have  converged  when  the  sublayers  totaled  14 
sublayers  (2^0  +  4Base  +  ^subgrade  4stifr)j  no  further  change  was  visible  as  the  number  of 
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sublayers  exceeded  the  standard.  The  remainder  of  the  analyses  were  performed  with  a 
standard  scheme  consisting  of  18  sublayers  (2ac  +  4Base  +  4subgrade  + 

The  peak  deflections  (Figure  7-8h)  increased  hy  16%  when  the  sublayers  were 
doubled  from  the  initial  setting,  but  only  increased  by  4%  with  further  increases.  The 
mean  (Figure  7-8c)  remained  nearly  constant  for  all  layering  combinations.  The  standard 
deviation  decreased  initially  by  7%,  but  remained  constant  for  combinations  with  more 
than  14  layers.  The  coefficient  of  skewness  (Figure  7-8d)  also  increased  initially,  but 
then  leveled  out,  indicating  that  the  basin  is  slightly  deeper  behind  the  mean  (relative  to 
the  direction  of  travel).  The  coefficient  of  kurtosis  (Figure  7-8e)  increased  by  10%  with 
the  first  increase  in  layers,  but  then  quickly  became  nearly  constant  when  more  than 
about  10  layers  were  used.  The  coefficient  of  kurtosis  indicates  that  the  shape  of  the 
deflection  basin  was  more  sharply  peaked  than  the  shape  of  a  normal  probability 
distribution. 

This  analysis  shows  that  the  accuracy  of  the  deflection  basin  is  not  improved  by 
increasing  the  number  of  layers  above  the  standard  scheme.  Based  on  this  observation 
and  the  significant  increase  in  computation  time  associated  with  more  layers,  the  standard 
layering  scheme  was  deemed  acceptable  and  used  for  the  remainder  of  the  analyses. 

7.2.2  Cutoff  Frequency 

As  described  in  Section  4.2.3,  the  loading  function  is  represented  by  a  bandwidth- 
limited  discrete  Fourier  series.  Increasing  the  cutoff  frequency  used  to  characterize  the 
load  improves  the  accuracy  of  the  representation,  but  also  increases  the  required 
computation  time. 

To  determine  the  best  cutoff  frequency  for  all  other  analyses,  six  specified  cutoff 
frequencies  were  examined,  each  with  a  frequency  increment  of  1 .95  Hz.  The  smallest 
cutoff  frequency  was  37.05  Hz  and  the  largest  was  154.05  Hz. 
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Figure  7-9  illustrates  the  calculated  basins  for  all  cutoff  frequencies.  The  lowest 
cutoff  frequency  (37.05  Hz)  resulted  in  a  slightly  deeper  and  wider  basin  than  higher 
cutoff  frequencies.  Virtually  no  change  was  recognizable  when  the  cutoff  frequency  was 
increased  above  76.05  Hz. 

The  peak  deflection  decreased  by  2%  when  the  cutoff  frequency  was  increased  from 
36  to  56  Hz,  but  remained  constant  thereafter.  The  mean  was  constant  for  all  cutoff 
frequencies.  The  standard  deviation  was  essentially  constant  at  1.17  until  the  cutoff 
frequency  exceeded  76  Hz,  and  then  jumped  by  8%  to  1.25  and  remained  constant 
thereafter.  The  coefficient  of  skewness  slowly  decreased  by  10%  as  the  cutoff  frequency 
increased  from  37  to  154  Hz.  The  coefficient  of  kurtosis  also  decreased  by  about  12% 
overall,  including  an  8%  drop  when  the  cutoff  frequency  increased  from  76.05  to  95.55 
Hz. 

This  analysis  showed  that  using  all  but  the  smallest  cutoff  frequency  produces 
approximately  the  same  deflection  basin.  Although  the  coefficient  of  skewness  and 
kurtosis  change  slightly  as  the  cutoff  frequency  was  increased,  unless  they  are  used  for 
backcalculation  purposes,  a  bandwidth  of  56  or  76  Hz  appears  to  be  acceptable.  For  this 
study,  76.05  Hz  was  selected  as  the  standard  cutoff  frequency  and  used  for  the  remainder 
of  the  calculations. 


7.2.3  Frequency  Increment 

The  frequency  increment,  hf,  used  in  a  discrete  Fourier  transformation  is 
determined  by  the  number  of  discrete  points,  N,  and  the  time  increment.  At,  between  each 
point  in  the  loading  function.  These  values  are  related  in  the  following  manner: 


Thus,  the  frequency  increment  can  be  increased  by  decreasing  the  time  increment  with 
the  total  number  of  points  held  constant  or  by  decreasing  the  total  number  of  points  with 
the  time  increment  held  constant.  Increasing  the  frequency  increment  decreases 
computation  time,  but  may  increase  error  in  the  calculated  deflection  basins. 


Peak  Deflection  [^m]  Deflection  [^m] 
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Figure  1-9.  Sensitivity  to  Cutoff  Frequency. 
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To  determine  the  best  frequency  increment,  three  combinations  of  Af  and  At  were 
used,  while  N  was  held  constant  at  1024. 

A/ =  0.00025  =>  A/ =  3.906 
At  =  0.0005  =>  A/  =  1.95 
Ar  =  0.001  =>  A/ =0.977 

Figure  7-10  presents  the  results  for  each  of  these  frequency  increments.  The  only 
visible  difference  occurred  when  the  frequency  increment  was  3.906  Hz.  For  this  case, 
the  basin  was  slightly  shallower  than  the  other  basins  at  a  distance  1  to  2  m  behind  the 
wheel.  Frequency  increments  of  1.95  and  0.977  Hz  produced  virtually  identical  basins. 

The  peak  deflection  for  each  frequency  increment  was  identical.  The  mean 
variation  showed  a  slight,  but  consistent  downward  trend;  the  mean  decreased  by  8% 
when  the  standard  frequency  increment  doubled  and  increased  by  8%  when  the  increment 
was  halved.  The  standard  deviation  increased  only  slightly  as  the  frequency  increment 
was  increased.  The  coefficient  of  skewness  was  largest  for  the  standard  increment,  and 
reduced  by  about  30%  for  the  other  two  increments.  The  coefficient  of  kurtosis  was 
lowest  for  the  standard  increment  and  increased  by  only  4%  and  1%  for  the  mininmim  and 
maximum  increments. 

Even  though  the  basins  did  not  change  visibly  for  these  different  frequency 
increments,  the  statistical  parameters  show  that  small  changes  did  occur.  The  minimum 
frequency  increment  should  produce  the  most  accurate  solution  because  the  loading  is 
most  accurately  represented.  However,  the  calculation  time  is  inversely  proportional  to 
the  frequency  increment,  i.e.  it  doubles  each  time  the  frequency  increment  is  halved. 
Therefore,  the  middle  increment  (1.95  Hz)  was  selected  as  the  standard  because  it 
appeared  to  produce  adequate  accuracy  in  a  reasonable  computation  time. 
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Figure  7-10.  Sensitivity  to  Frequency  Increment. 


7.2.4  Timeback 
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The  variable  timeback  designates  how  far  back  in  time  (and  space)  the  load  pulses 
are  considered  in  deflection  calculations.  An  appropriate  value  for  timeback  could  be 
selected  by  closely  examining  the  deflection  time  histories  for  load  pulses  at  various 
distances  and  determining  when  the  response  has  ceased.  Alternately,  the  value  of 
timeback  could  be  adjusted  to  determine  the  actual  influence  on  rolling- wheel  deflection 
basins.  Like  the  other  parameters  examined  in  this  section,  timeback  has  a  strong 
influence  on  the  required  computation  time  —  the  number  of  load  pulses  for  which  a 
response  must  be  calculated  is  directly  proportional  to  timeback. 

Deflection  basins  using  four  different  values  of  timeback  were  calculated,  ranging 
from  0.10  to  0.25  seconds.  For  the  standard  vehicle  speed  of  50  km/hr  =  13.89  m/s  (31 
mph),  these  timeback  values  correspond  to  distances  of  1 .39  to  3.47  m  (4.6  to  1 1 .4  ft). 

Figure  7-11  illustrates  the  calculated  basins  for  all  timeback  values.  The  only 
discernible  change  in  the  basins  appeared  at  distances  between  1  and  3  m  (3.3  and  9.8  ft) 
behind  the  wheel.  The  basin  became  very  slightly  deeper  in  this  zone  as  timeback  was 
increased  from  0.15  to  0.20  sec.  No  visible  change  was  observed  when  timeback  wus 
extended  to  0.25  sec. 

The  peak  deflection  was  constant  for  all  values  of  timeback.  The  mean  increased 
very  slightly  as  timeback  increased,  and  the  standard  deviation  decreased  by  3% 
when  timeback  was  doubled  from  0.10  to  0.20.  The  coefficient  of  skewness  was  largest 
when  timebackwas  equal  to  0.15,  then  dropped  by  12%  as  timeback  increased  to  0.25. 
The  coefficient  of  kurtosis  increased  by  only  1.5%  when  timeback  was  doubled  from  0.1 
to  0.2  see. 

This  analysis  shows  that  the  deflection  basin  did  not  change  visibly  as  timeback  is 
increased  beyond  0.2  sec.  The  standard  deviation  and  mean  did  change  slightly  as 
timeback  was  increased.  Additionally,  the  computation  time  increased  by  about  50% 
when  timeback  was  increased  from  0.20  to  0.25  msec.  Therefore,  the  standard  value  for 
timeback  was  set  at  0.20  seconds  and  used  for  the  remainder  of  the  analyses. 


7.2.5  Avetime 
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The  parameter  avetime  represents  the  time  step  between  each  discrete  load  pulse. 

Because  the  response  to  each  pulse  must  be  computed  separately,  the  number  of  pulses, 

N,  considered  within  timeback  has  a  significant  influence  on  the  required  computation 

time.  Avetime  and  timeback  are  related  in  the  following  manner. 

timeback 

N  = - 

avetime 

To  determine  the  best  time  increment  between  load  pulses,  four  different  values  of 
avetime  were  analyzed,  ranging  from  0.5  to  4.0  msec.  For  the  standard  vehicle  speed  of 
50  km/hr  =  13.89  m/s  (31  mph),  these  avetime  values  correspond  to  distances  of  7  to  56 
mm  (0.3  to  2.2  in). 

Figure  7-12  illustrates  the  calculated  basins  for  all  values  of  avetime.  When  avetime 
was  equal  to  4  msec,  the  basin  was  shifted  slightly  forward  of  the  other  three  basins, 
however,  the  peak  deflection  was  nearly  identical  to  the  others. 

The  peak  deflection  for  the  three  smallest  values  of  avetime  was  nearly  identical, 
and  when  avetime  was  increased  from  1  to  4  msec,  the  peak  decreased  by  only  0.7%.  The 
mean  was  constant  for  the  two  smallest  values  of  avetime,  but  decreased  by  8%  as 
avetime  was  increased  to  2  and  4  msec.  The  coefficient  of  skewness  showed  a  slightly 
downward  trend,  decreasing  a  total  of  20%  from  the  minimum  to  maximimi  avetime.  The 
coefficient  of  kurtosis  was  nearly  constant  for  all  avetime. 

In  general,  the  basin  shapes  and  the  statistical  parameters  converged  to  constant 
values  as  avetime  was  decreased.  The  computation  time  tripled  when  avetime  was 
changed  from  1.0  to  0.5  msec,  with  no  significant  improvement  in  accuracy.  Therefore, 
the  standard  avetime  value  was  set  at  1  msec  and  used  for  the  remainder  of  the  analyses. 
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Figure  7-12.  Sensitivity  to  Avetime. 
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7.3  Material  Variables 

The  material  variables  describe  the  physical  characteristics  of  the  pavement 
structure.  They  are  the  properties  that  pavement  tests  are  designed  to  determine  because 
they  can  provide  an  indication  of  the  strength  and  life  expectancy  of  a  pavement  system. 
Current  backcalculation  schemes  can  only  determine  one  or  two  types  of  these  variables, 
so  the  others  must  be  assumed  or  measured.  Typically,  the  density,  Poisson’s  ratio,  and 
thickness  of  each  layer  are  predetermined  and  the  modulus  is  backcalculated.  This 
section  describes  the  deflection  basin  sensitivity  to  each  of  these  material  variables. 

7.3.1  Density 

The  density  of  each  pavement  layer  is  most  strongly  related  to  the  type  and 
gradation  of  the  material  in  each  layer.  AC  layers  have  the  highest  densities  because  they 
are  primarily  composed  of  crushed  stone  aggregate  with  the  voids  between  aggregate 
particles  filled  with  asphalt  cement.  Subgrade  layers  generally  have  a  lower  density 
because  they  are  often  composed  of  fine-grained  soil  with  very  little  aggregate.  The 
gravel  mixtures  used  in  base  course  materials  usually  produce  a  density  somewhere 
between  the  AC  and  subgrade. 

Three  density  schemes  were  selected  for  testing  and  applied  to  all  material  layers: 
10%  less  than  standard,  standard,  and  10%  greater  than  standard.  Material  layer  densities 
can  usually  be  estimated  within  a  10  or  20%  error,  so  the  limited  range  of  densities  was 
reasonable. 

Figure  7-13  illustrates  the  calculated  basins  for  all  density  combinations.  There  was 
no  visible  change  in  the  deflection  basins  for  either  offset  combination  of  densities. 

The  peak  deflection  for  both  offset  density  combinations  was  0.4%  less  than  that  of 
the  standard  densities.  The  mean  and  standard  deviation  were  essentially  constant  for  all 
density  combinations.  The  coefficient  of  skewness  increased  by  3%  from  the  standard  for 
both  of  the  other  density  combinations.  The  coefficient  of  kurtosis  decreased  by  only 
0.4%  from  the  lowest  to  the  highest  density  combination. 
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This  analysis  shows  that  the  deflection  basin  is  extremely  insensitive  to  the 
density  of  the  material  layers.  Therefore,  even  if  the  densities  are  slightly  in  error,  the 
predicted  deflection  basins  will  still  be  accurate. 

7.3.2  Poisson’s  Ratio 

Poisson’s  ratio  relates  the  amount  of  lateral  strain  to  the  axial  strain  resulting  from 
loading  in  the  axial  direction.  Poisson’s  ratio,  v,  also  relates  Yoxmg’s  modulus,  E,  to 
shear  modulus,  G,  in  an  isotropic  linear  elastic  material  in  the  following  manner. 

G=  , 

2(1 +v) 

Five  combinations  of  Poisson’s  ratio  in  the  layers  were  tested,  as  listed  below. 


Laver 
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0.325 

0.35 

0.35 

0.35 

0.35 
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0.35 

0.35 

0.40 

0.40 

0.40 

Subgrade 

0.35 

0.35 

0.40 

0.45 

0.45 

Stiff  Layer 

0.35 

0.35 

0.40 

0.35 

0.45 

The  Poisson’s  ratio  for  each  layer  is  usually  estimated  for  backcalculation,  but  there  is 
considerable  range  in  the  values  that  are  commonly  used.  The  values  listed  above  were 
selected  to  represent  typical  values  used  for  pavement  deflection  backcalculation. 
Combination  4  was  selected  as  the  standard  setting. 

Figure  7-14  illustrates  the  calculated  basins  for  all  Poisson’s  ratio  combinations. 
There  was  very  little  visible  change  in  the  basins;  however,  close  examination  reveals 
that  the  first  three  combinations  resulted  in  slightly  deeper  and  wider  deflection  basins. 
There  was  no  visible  difference  between  Combinations  4  and  5. 


Peak 


Figure  7-14.  Sensitivity  to  Poisson’s  Ratio. 
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The  peak  deflection  for  Combination  1  was  4%  greater  than  that  of  the  standard, 
but  all  other  combinations  were  within  2%  of  the  standard.  The  mean  was  nearly 
constant  for  all  combinations  and  Combination  2  produced  the  largest  difference  in 
standard  deviation  —  6%  greater  than  the  standard.  The  coefficient  of  skewness  was  12% 
less  than  the  standard  for  Combinations  1  and  2,  but  closer  to  standard  for  Combinations 
3  and  5,  at  7%  less  and  1%  greater,  respectively.  The  coefficient  of  kurtosis  was  more 
than  12%  less  than  standard  for  Combinations  1  and  2,  and  again.  Combinations  3  and  5 
were  closer  to  standard. 

This  analysis  shows  that  the  deflection  basin  is  fairly  insensitive  to  the  Poisson’s 
ratio  of  the  each  layer.  The  tested  range  covered  the  extremes  of  assumed  values,  yet 
none  of  the  parameters  was  more  than  14%  different  than  the  standard.  This  indicates 
that  Poisson’s  ratios  should  be  carefully  estimated  but  large  errors  will  not  result  if  the 
selected  values  are  slightly  in  error. 

7.3.3  Modulus 

The  modulus  (or  stiffiiess)  of  a  material  relates  applied  stress  to  resulting  strain.  In 
pavement  testing  and  analysis,  the  modulus  of  each  layer  strongly  influences  the 
deflections  produced  by  a  wheel  load.  Large  deflections  result  from  large  strains  which 
ultimately  cause  the  pavement  materials  to  weaken  and  deteriorate.  This  section 
describes  the  deflection  basin  sensitivity  to  the  modulus  of  each  pavement  layer. 

7.3.3. 1  AC  Modulus 

The  modulus  of  the  AC  layer  was  assigned  eight  values  between  250  and  2500  ksi 
(1724  and  17238  MPa),  with  the  standard  AC  modulus  set  at  562.8  ksi  (620.6  MPa).  AC 
moduli  can  exceed  1000  ksi  (6895  MPa)  at  asphalt  temperatures  below  10°  C  (50°  F).  At 
the  other  end  of  the  scale,  an  AC  modulus  of  250  ksi  can  correspond  to  a  distressed  or  a 
very  warm  pavement  (>27°  C  (80°  F)). 
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Figure  7-15  illustrates  the  deflection  behavior  calculated  for  each  AC  modulus. 

The  majority  of  the  changes  in  the  basin  occurred  within  0.8  m  (2.6  ft)  of  the  peak 
deflection,  and  deflections  farther  from  the  peak  changed  very  little  over  the  entire  range 
of  AC  moduli.  The  basin  became  very  pointed  when  the  AC  modulus  was  small,  and 
much  more  rounded  when  the  AC  was  extremely  large.  Close  observation  of  the  area 
about  1  m  ahead  and  behind  the  peak  showed  that  the  deflections  increased  in  this  zone  as 
the  AC  stiffness  increased,  resulting  in  a  wider  deflection  basin. 

The  peak  deflection  decreased  by  20%  when  the  modulus  was  doubled  from  the 
initial  value,  and  subsequent  doublings  resulted  in  22%  and  23%  reductions  from  the 
previous  value.  The  mean  was  nearly  constant  over  the  entire  range  of  moduli,  while  the 
standard  deviation  steadily  increased  as  the  AC  modulus  increased,  registering  an  1 8% 
increase  from  the  smallest  to  largest  values.  The  coefficient  of  skewness  slowly 
decreased  by  about  20%  over  the  range  of  moduli,  indicating  that  the  deflection  basin 
became  more  symmetric  with  increasing  AC  modulus.  The  coefficient  of  kurtosis  rapidly 
decreased  as  the  modulus  increased,  which  should  be  expected  from  observing  the 
relative  peakedness  of  the  basins. 

1332  Base  Modulus 

The  modulus  of  the  base  layer  was  varied  between  5  and  50  ksi  (35  and  345  MPa), 
with  the  standard  layer  modulus  set  at  14.8  ksi  (102  MPa).  A  base  modulus  of  5  ksi 
corresponds  to  a  relatively  soft  fine-grained  soil  and  a  modulus  of  50  ksi  may  be  achieved 
with  a  well-graded,  well-compacted  crushed  stone  mixture.  Most  base  course  materials 
have  a  modulus  between  12  and  30  ksi  (83  and  207  MPa). 

Figure  7-16  illustrates  the  deflection  behavior  calculated  for  each  base  modulus 
value.  The  base  course  modulus  influenced  deflections  within  about  1.2  m  (3.9  ft)  of  the 
peak  deflection.  Almost  no  change  in  the  basins  was  visible  outside  this  zone. 


Ill 

The  peak  deflection  increased  by  20%  when  the  base  modulus  was  reduced  from 
the  standard  by  65%.  The  peak  deflection  decreased  by  28%  when  the  base  modulus  was 
increased  from  the  standard  to  the  maximum.  The  mean  value  was  essentially  constant 
for  all  beise  moduli,  while  the  standard  deviation  increased  by  20%  from  the  minimum  to 
maximum  modulus  tested.  The  coefficient  of  skewness  slowly  decreased  by  15%  over 
the  range  of  moduli,  indicating  that  the  deflection  basin  becomes  more  symmetric  as  the 
modulus  increases.  The  coefficient  of  kurtosis  increased  by  18%  from  the  standard  to  the 
minimum,  but  only  decreased  by  19%  from  the  standard  to  the  maximum  modulus. 

7.3. 3. 3  Subgrade  Modulus 

The  modulus  of  the  subgrade  layer  was  varied  between  3  and  25  ksi  (21  and  173 
MPa),  with  the  standard  layer  modulus  set  at  10.2  ksi  (70.2  MPa).  A  subgrade  modulus 
of  3  ksi  corresponds  to  a  very  soft  clay  or  very  loose  sand,  and  a  modulus  of  25  ksi  to  a 
very  dense  well-graded  soil.  The  range  of  subgrade  moduli  typically  encountered  is 
highly  dependent  on  the  local  soil  conditions.  Desirable  subgrades  for  highway 
construction  generally  have  moduli  greater  than  about  12  ksi  (83  MPa). 

Figure  7-17  illustrates  the  deflection  behavior  calculated  for  each  value  of  subgrade 
modulus.  The  subgrade  modulus  influences  deflections  over  a  wide  area  of  the 
pavement,  stretching  nearly  2  m  (6.6  ft)  in  either  direction  from  the  wheel.  The  general 
shape  of  the  basin  stayed  fairly  consistent  in  the  upper  range  of  moduli.  Heaving 
(positive  deflections)  was  predicted  about  3  m  (9.8  ft)  behind  the  wheel  for  very  weak 
moduli  of  3  ksi  (21  MPa)  or  less. 

The  peak  deflection  rapidly  decreased  as  the  modulus  increased  from  the  minimum 
to  the  standard.  Above  the  standard,  the  peak  deflections  decreased  nearly 
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Figure  7-17.  Sensitivity  to  Subgrade  Modulus. 
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linearly,  declining  about  4  |xm  for  every  1  ksi  increase  in  modulus.  The  mean 
increased  slightly  over  the  3  to  10  ksi  range,  but  then  became  constant  as  the  modulus 
increased.  The  standard  deviation  increased  by  20  %  when  the  modulus  was  doubled 
from  5  to  10  ksi,  but  then  increased  less  rapidly  from  the  standard  to  the  Tnayimnm 
modulus.  The  coefficient  of  skewness  sharply  dropped  by  nearly  40%  when  the  subgrade 
modulus  was  increased  from  3  to  5  ksi;  thereafter,  it  decreased  slowly  over  the  middle 
range  of  moduli.  The  coefficient  of  kurtosis  posted  a  15%  gain  when  the  modulus 
increases  from  3  to  5,  but  steadily  decreased  as  the  modulus  rose  to  25  ksi. 

7.3.3 .4  Stiff  Layer  Modulus 

The  modulus  of  the  stiff  layer  was  varied  between  5  and  1000  ksi,  with  the  standard 
modulus  set  at  40  ksi.  A  modulus  of  5  ksi  corresponds  to  a  soft  clay  or  very  loose  sand 
deposit  —  not  exactly  a  “stiff’  layer,  but  thick  deposits  of  these  materials  do  exist. 

Moduli  approaching  1000  ksi  may  represent  bedrock  or  a  frozen  groundwater  layer. 

Figure  7-18  illustrates  the  deflection  behavior  calculated  for  each  stiff  layer 
modulus.  The  stiff  layer  modulus  had  a  wider  influence  on  the  deflection  basin  than  any 
other  layer  modulus.  Visible  deflections  extended  beyond  5  meters  from  the  wheel  when 
the  modulus  was  equal  to  5  and  10  ksi.  As  the  stiff  layer  modulus  was  increased  above 
25  ksi,  the  basins  began  to  converge  to  a  consistent  shape.  Heave  was  predicted  about 
1 .5m  in  front  of  the  wheel  and  2m  behind  the  wheel  when  the  modulus  reached  1000  ksi. 

The  peak  deflection  initially  decreased  rapidly  by  14%  when  the  modulus  was 
doubled  from  5  to  10  ksi;  above  25  ksi,  the  decrease  was  nearly  linear.  The  mean 
noticeably  increased  as  the  modulus  increased  from  5  to  40  ksi,  and  then  leveled  out 
thereafter.  The  standard  deviation  was  very  sensitive  to  the  stiff  layer  modulus, 
decreasing  by  65%  over  the  range  from  5  to  75  ksi.  The  coefficient  of  skewness  steadily 
increased  over  the  range  of  moduli  examined.  The  coefficient  of  kurtosis  rapidly 
increased  as  the  modulus  increased,  gaining  245%  from  the  minimum  to  the  maximum 
modulus  examined. 
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7.33.5  Discussion 

The  deflection  basin  produced  by  a  rolling  wheel  is  very  sensitive  to  the  modulus  of 
each  material  layer.  As  expected,  the  modulus  of  the  deeper  material  layers  influenced  a 
wider  zone  of  the  deflection  basin.  Very  soft  subgrade  and  stiff  layers  produce  such  a 
wide  deflection  basin  that  they  may  be  very  difficult  to  measure  because  the  imdeflected 
reference  measurement  could  easily  be  inside  the  deflection  basin.  Overall,  the  peak 
deflection  and  the  standard  deviation  are  the  most  consistently  sensitive  statistical 
parameters  for  each  layer  modulus. 

7.3.4  Thickness 

The  thickness  of  each  pavement  layer  is  a  key  design  decision  when  constructing 
new  pavements.  For  pavement  testing  and  analysis,  the  thickness  is  an  important  input 
for  parameter  backcalculation  of  layer  moduli.  Typically,  the  layer  thickness  is 
determined  from  drilled  cores,  construction  specifications,  or  occasionally  field  tests. 

The  sensitivity  of  the  dynamic  deflection  basin  to  layer  thickness  may  be  such  that  the 
layer  thicknesses  can  be  backcalculated  with  the  layer  moduli.  This  section  investigates 
the  deflection  basin  sensitivity  to  the  thickness  of  each  layer. 

7.3.4. 1  AC  Thickness 

The  thickness  of  the  AC  layer  was  varied  between  50  and  600  mm  (2  and  24  in), 
with  the  standard  layer  thickness  at  137  mm  (5.4  in).  AC  layers  as  thin  as  50  mm  may  be 
used  for  a  temporary  road,  and  600  mm  of  asphalt  could  be  specified  for  a  cargo  aircraft 
runway.  Many  highway  AC  layers  are  between  150  and  300  mm  (5.9  and  1 1.8  in)  thick. 

Figure  7-19  illustrates  the  computed  deflection  behavior  for  each  AC  thickness. 

The  largest  changes  in  the  basin  occurred  within  0.8  m  (2.6  ft)  of  the  peak  deflection, 
where  the  deflections  decreased  as  the  AC  thickness  increased.  However,  the  deflections 
outside  this  zone  actually  increased  as  the  AC  thickness  increased.  At  large  thicknesses, 
the  basin  began  to  show  a  distinctive  roxmded  peak  within  0.2  m  (0.7  ft)  of  the  peak 
deflection;  beyond  the  wheel  contact  area,  the  basin  became  very  broad. 
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Figure  7-19.  Sensitivity  to  AC  Thickness. 
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The  peak  deflection  decreased  rapidly  as  the  thickness  increased  from  50  to  200 
mm,  but  then  decreased  by  64%  when  the  thickness  was  tripled  from  200  to  600  mm. 

The  mean  decreased  slightly  from  the  minimum  to  maximum  thickness,  while  the 
standard  deviation  increased  nearly  linearly  with  increasing  AC  thickness.  The 
coefficient  of  skewness  decreased  by  about  0.06  for  every  100  mm  increase  in  thickness. 
In  a  manner  consistent  with  the  shapes  of  the  basins,  the  coefficient  of  kurtosis  rapidly 
decreased  as  the  thickness  increased. 

7.3. 4.2  Base  Thickness 

The  thickness  of  the  base  layer  was  varied  between  200  and  800  mm  (8  and  31.5 
in),  with  the  standard  layer  thickness  at  330  mm  (13  in).  A  200  mm  thickness  can  be 
used  when  the  subgrade  is  fairly  competent,  while  a  thickness  of  800  mm  might  be 
required  when  the  subgrade  is  weak  or  soft.  Typical  base  thicknesses  range  from  1 50  to 
600  mm  (6  to  24  in). 

Figure  7-20  illustrates  the  deflection  behavior  calculated  for  each  base  thickness. 
The  only  visible  changes  in  the  deflection  basins  occurred  in  the  zones  1  to  2  m  (3.3  to 
6.6  ft)  ahead  and  behind  the  peak  deflections.  Very  small  changes  were  also  observed  at 
the  peak  of  the  basin. 

The  peak  deflection  was  nearly  constant,  only  decreasing  3%  over  the  entire  range 
of  thicknesses  studied.  The  mean  was  also  nearly  constant;  the  standard  deviation  was 
constant  for  thicknesses  between  200  and  600  mm,  but  then  jumped  by  9%  when  the 
thickness  was  increased  from  600  to  800  mm.  The  coefficient  of  skewness  was  nearly 
constant,  with  a  slight  dip  observed  for  the  400  mm  thickness.  The  coefficient  of  kurtosis 
decreased  nearly  linearly  by  0.056  for  every  100  mm  increase  in  thickness. 
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Figure  7-20.  Sensitivity  to  Base  Thickness. 
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7.3. 4.3  Subgrade  Thickness 

The  thickness  of  the  subgrade  layer  was  varied  between  0.5  and  2.0  m  (20  and  79 
in),  with  the  standard  layer  thickness  at  1.17  m  (46  in).  The  subgrade  thickness  of  an 
actual  pavement  system  depends  on  the  local  site  conditions.  The  bottom  of  the  subgrade 
is  defined  by  the  top  of  the  stiff  layer;  therefore,  the  thickness  is  entirely  determined  by 
the  depth  to  bedrock  or  the  ground  water  table. 

Figure  7-21  illustrates  the  deflection  behavior  calculated  for  each  subgrade 
thickness.  As  the  subgrade  thickness  increased,  the  deflection  basin  became  wider  and 
deeper,  because  the  underlying  stiff  layer  was  deeper  below  the  pavement  surface. 

The  peak  deflection  decreased  by  10%  when  the  thickness  was  halved  from  1 .0  to 
0.5  m,  and  increased  by  10%  when  the  thickness  was  doubled  from  1 .0  to  2.0  m.  The 
mean  was  also  nearly  constant  over  the  range  of  thicknesses;  the  standard  deviation 
decreased  by  10%  when  the  thickness  was  halved  from  1.0  to  0.5  m,  but  decreased  by 
only  7.5%  when  the  thickness  doubled  from  1 .0  to  2.0  m.  The  coefficient  of  skewness 
was  largest  for  the  standard  thickness,  and  decreased  by  32%  when  the  thickness  was 
increased  to  2.0  m.  The  coefficient  of  kurtosis  showed  the  same  trend  as  the  skewness, 
with  the  standard  thickness  exhibiting  the  largest  kurtosis. 

7.3. 4.4  Stiff  Layer  Thickness 

The  thickness  of  the  stiff  layer  was  varied  between  5  and  40  m  (16.4  and  131  ft), 
with  the  standard  layer  thickness  at  10  m  (33  ft).  Because  the  material  below  the  stiff 
layer  is  modeled  as  a  halfspace  with  the  same  properties  as  the  stiff  layer,  the  thickness  of 
the  subgrade  only  needs  to  be  large  enough  to  ensure  an  accurate  solution. 

Figure  7-22  illustrates  the  deflection  behavior  calculated  for  each  stiff  layer 
thickness.  Only  the  basin  calculated  for  a  stiff  layer  5  m  thick  was  visibly  different  than 
the  other  basins. 
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Figure  7-21.  Sensitivity  to  Subgrade  Thickness. 
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The  peak  deflection  for  the  5  m  layer  was  2.5%  less  than  the  deflection  for  the 
standard  (and  all  other)  thickness.  The  mean  was  constant  for  all  cases;  the  standard 
deviation  was  30%  less  than  that  of  the  standard  when  the  stiff  layer  thickness  was  5  m, 
and  only  increased  3%  above  standard  when  the  stiff  layer  thickness  was  increased  to  40 
m.  The  coefficient  of  skewness  decreased  by  13%  when  the  thickness  was  increased 
from  5  to  20  m,  but  was  constant  for  thicker  layers.  The  coefficient  of  kurtosis  was 
largest  for  the  standard  thickness,  and  decreased  by  5%  when  the  thickness  was  increased 
to  40  m. 

7.3. 4. 5  Discussion 

The  deflection  basin  produced  by  a  rolling  wheel  is  very  sensitive  to  the  thickness 
of  the  AC  layer,  and  somewhat  sensitive  to  the  thickness  of  the  subgrade  layer.  However, 
the  base  thickness  had  very  little  influence  on  the  deflection  basin.  The  standard 
deviation  and  coefficient  of  kurtosis  were  very  sensitive  to  the  AC  thickness,  but  not  to 
the  thicknesses  of  the  other  layers.  The  standard  stiff  layer  thickness  of  10  m  was  shown 
to  be  an  acceptable  choice  because  the  basin  shape  did  not  visibly  change  as  the  thickness 
was  increased. 

7.3.5  Combinations  of  Thickness  and  Modulus 

Most  pavement  backcalculation  programs  cannot  determine  both  the  thickness  and 
modulus  of  each  layer.  Generally  the  thickness  is  measured  or  assumed,  and  the  modulus 
is  backcalculated.  This  section  examines  the  sensitivity  of  the  deflection  basin  to  both 
modulus  and  thickness. 

The  static  backcalculation  program  EVERCALC  4.0  was  used  to  determine  the 
moduli  of  pavement  layers  for  given  input  thicknesses.  By  specifying  the  same  static 
deflection  basin  while  adjusting  the  thickness  of  only  one  layer,  the  modulus  of  all  of  the 
layers  was  backcalculated.  With  this  scheme,  a  stiff,  thin  layer  and  a  soft,  thick  layer 
would  produce  identical  static  deflection  basins.  In  all  cases,  the  stiff  layer  was  held 
constant  at  1.64  m  (64.4  in)  below  the  pavement  surface. 
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7.3.5. 1  AC  Thickness/Modulus 

Five  combinations  of  AC  thickness  and  modulus,  listed  in  Table  7-1,  were 
examined.  The  modulus  of  the  base  and  subgrade  layers  was  adjusted  for  each  case  to 
produce  identical  static  deflection  basins. 


Table  7-1.  AC  Thickness/Modulus  Combinations. 


-50% 

-25% 

Standard 

10% 

20% 

AC  Thickness  [mm] 

69 

103 

137 

171 

206 

AC  Modulus  [MPa] 

15.1 

4.7 

2.2 

1.3 

0.9 

Figure  7-23  illustrates  the  deflection  behavior  calculated  for  each  AC 
thickness/modulus  combination.  The  only  visible  changes  in  the  deflection  basin  were 
observed  tvithin  about  0.2  m  (0.7  ft)  of  the  peak  deflection.  In  this  zone,  the  peak 
deflection  increased  with  increasing  AC  thickness. 

The  peak  deflection  steadily  increased  as  the  AC  thickness  increased,  gaining  17% 
from  the  minimum  to  the  maximum  thickness.  The  mean  and  standard  deviation  did  not 
change  appreciable  over  the  range  of  thicknesses.  The  coefficient  of  skewness  increased 
slightly  with  AC  thickness,  gaining  6%  from  the  minimum  to  maximum  thickness.  The 
coefficient  of  kurtosis  was  nearly  constant  for  all  but  the  thinnest  AC,  which  was  2.7% 
less  than  the  standard. 

7.3. 5.1  Base  Thickness/Modulus 

Five  combinations  of  base  thickness  and  modulus,  listed  in  Table  7-2,  were 
examined.  The  modulus  of  the  AC  and  subgrade  layers  was  adjusted  for  each  case  so  that 
the  static  deflection  basin  remained  constant. 


Peak  Deflection  [  ^m]  Deflection  [^m] 


125 


Distance  [m] 

-5  -4  -3  -2  -1  0  1  2  3  4  5 


O 


0.6 


0.4 


0.2 


0.0 


8  - 

- - - ♦ - e - 0 

o - e - ♦ - ^ ^ 

6  - 

O 

4  - 

2  - 

- 1 - 1 - ^ ^ - 

0  - 

- i - i - 1 - ^ - 

50  100  150  200 

AC  TTiickness  Combination  [mm] 


250 


50  100  150  200 

AC  ITiickness  Combination  [mm] 


250 


Figure  7-23.  Sensitivity  to  AC  Thickness/Modulus. 
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Table  7-2.  Base  Thickness/Modulus  Combinations. 


-50% 

-25% 

Standard 

10% 

20% 

Base  Thickness  [mm] 

165 

247 

330 

413 

495 

Base  Modulus  [MPa] 

0.343 

0.162 

0.106 

0.101 

0.092 

Figure  7-24  illustrates  the  deflection  behavior  calculated  for  each  base 
thickness/modulus  combination.  In  this  case,  the  deflections  only  changed  significantly 
within  0.1  m  (3.9  in)  of  the  peak.  As  the  base  thickness  increased,  the  peak  deflection 
decreased,  most  significantly  when  the  base  was  50%  thinner  than  the  standard. 

The  peak  deflection  was  non-linearly  related  to  the  base  thickness.  The  peak 
deflection  for  the  thinnest  base  was  7%  greater  than  the  standard,  but  the  peak  deflection 
for  the  thickest  base  was  only  2%  less  than  the  standard.  The  mean  and  standard 
deviations  were  nearly  constant  for  all  thickness/modulus  combinations.  The  coefficient 
of  skewness  decreased  slightly  as  the  base  thickness  increased,  declining  only  5%  overall. 
The  coefficient  of  kurtosis  also  decreased  over  the  range  of  base  thickness/modulus 
combinations,  declining  only  4%  overall. 

7.3. 5. 2  Discussion 

The  changes  in  deflection  basins  for  both  analyses  were  surprisingly  small, 
considering  the  wide  range  of  thicknesses  used.  The  AC  and  base  layer  showed  opposite 
trends  for  the  peak  deflection,  skewness,  and  kurtosis  change  with  thickness.  For  the  AC 
layer,  all  three  parameters  increased  as  the  AC  thickness  increased.  For  the  base  layer,  all 
three  parameters  decreased  as  the  base  thickness  increased.  However,  the  absolute 
changes  in  basins  were  so  small  that  it  appears  that  backcalculating  both  the  modulus  and 
thickness  may  be  difficult.  The  results  suggest  that  any  attempt  at  simultaneous 
backcalculation  of  modulus  and  layer  thickness  will  require  accurate  definition  of  the 
deepest  part  of  the  deflection  basin. 
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7.4  Loading  Variables 

The  loading  variables  are  a  function  of  the  vehicle  speed,  weight,  and  tire  pressure. 
The  deflection  basin  sensitivity  to  these  parameters  may  help  determine  the  operating 
conditions  that  would  provide  the  most  useful  information  about  the  pavement. 
Additionally,  the  optimum  location  for  deflection  sensors  might  be  determined  from  by 
examining  these  parameters. 

7.4.1  Wheel  T.oad 

The  wheel  load  was  varied  between  3,000  and  1 1,000  lb  (13,344  and  48,928  N), 
with  the  standard  layer  thickness  at  5,540  lb  (24,642  N).  The  tire  pressure  was  held 
constant  at  90.0  psi  (620.6  kPa);  therefore,  the  loaded  area  was  adjusted  so  that  constant 
stress  was  applied  for  all  loads.  A  3,000  lb  load  corresponds  to  a  lightly  loaded  panel 
truck,  and  an  1 1,000  lb  load  represents  an  overloaded  wheel  (9,000  lb  is  the  limit  in  most 
states). 

Figure  7-25  illustrates  the  computed  deflection  behavior  for  each  wheel  load.  As 
the  load  increased,  the  basin  became  consistently  deeper  and  wider. 

The  peak  deflection  linearly  gained  about  50  pm  with  each  1,000  lb  increase  in 
wheel  load.  The  mean  was  constant  for  all  loads,  and  the  standard  deviation  was  nearly 
constant  as  well,  gaining  only  0.9%  from  the  minimum  to  maximum  load.  The 
coefficient  of  skewness  was  also  nearly  constant,  and  the  coefficient  of  kurtosis  showed  a 
slightly  decreasing  trend  as  the  load  increased. 

The  deflection  basin  behavior  as  the  load  increased  was  exactly  as  would  be 
expected  for  a  linear-elastic  analysis.  When  the  load  is  doubled,  the  deflections  at  all 
locations  double.  This  is  confirmed  by  consistency  of  all  the  statistical  parameters. 


7.4.2  Tire  Pressure 
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Four  tire  pressures  were  examined,  ranging  from  30  to  120  psi  (207  to  827  kPa). 

The  standard  value  of  90.0  psi  (620.6  kPa)  is  a  typical  pressure  for  heavy  trucks.  Lower 
truck  tire  pressures  are  occasionally  used  when  heavy  loads  are  carried  across  weak 
pavement  structures,  such  as  temporary  gravel  roads.  Higher  tire  pressures  are  used 
because  they  improve  gas  mileage  at  high  speeds.  In  each  case  examined  in  this  study,  a 
constant  wheel  load  was  applied  and  the  load  radius  was  adjusted  accordingly. 

Figure  7-26  illustrates  the  deflection  behavior  calculated  for  each  tire  pressure.  As 
the  load  increased,  the  basin  became  slightly  deeper,  but  only  the  area  within  about  0.2  m 
(7.9  in)  of  the  peak  deflection  changed.  This  zone  of  change  may  reflect  the  size  of  the 
load  radius.  The  radius  of  the  loaded  area  was  0.195  m  (7.7  in)  at  30  psi  (207  kPa),  but 
reduced  to  0.097  m  (3.8  in)  at  120  psi  (827  kPa).  These  radii  correspond  to  the  areas  of 
the  basins  that  actually  change  with  the  tire  pressure. 

The  peak  deflection  decreased  by  6%  when  the  tire  pressure  was  doubled  from  30  to 
60  psi,  but  only  increased  by  3%  when  the  pressure  was  doubled  again.  The  standard 
deviation  decreased  slightly  as  the  pressure  increased,  but  the  mean  remained  essentially 
constant.  The  skewness  also  was  constant  for  all  the  tire  pressures.  The  coefficient  of 
kurtosis  very  slightly  increased  as  the  tire  pressure  was  increased,  reflecting  the  smaller 
contact  area  at  higher  tire  pressures. 

The  deflection  basin  was  not  very  sensitive  to  tire  pressure.  However,  this  analysis 
does  show  that  higher  pressures  produce  larger  deflections  directly  under  the  loaded  area, 
supporting  the  use  of  low  tire  pressures  on  weak  roads. 

7.4.3  Vehicle  Speed 

The  wheel  speed  was  varied  between  10  and  90  km/hr  (6  and  56  mph).  The  speed 
was  adjusted  by  spacing  the  same  load  pulses  further  apart.  In  this  manner,  the  same 
loading  was  applied  over  time  for  all  speeds,  but  the  loading  was  significantly  changed 


over  space. 


132 


Figvire  7-27  illustrates  the  computed  deflection  behavior  for  each  wheel  speed. 

As  the  velocity  increased,  the  peak  deflection  steadily  reduced  and  the  location  of  the 
peak  moved  steadily  farther  behind  the  wheel.  At  the  slowest  speed,  the  basin  was  nearly 
symmetrical  underneath  the  wheel. 

The  peak  deflection  linearly  decreased  by  1 1  pm  for  each  10  km/hr  increase  in 
speed.  The  mean  was  directly  proportional  to  the  speed,  increasing  (moving  farther 
behind  the  wheel)  100  mm  for  each  10  km/hr  increase.  The  standard  deviation  increased 
slightly  with  speed,  showing  a  larger  increase  per  interval  at  higher  speeds.  The 
coefficient  of  skewness  increased  in  a  very  non-linear  fashion.  Initially,  the  coefficient  of 
skewness  was  only  0.08,  indicating  the  basin  was  nearly  symmetric  about  the  mean.  As 
the  speed  increased,  the  coefficient  of  skewness  increased,  but  less  with  each  increment, 
eventually  leveling  off  at  about  0.35.  The  coefficient  of  kurtosis  steadily  decreased  about 
0.2  for  each  incremental  increase  in  speed. 

For  this  variable,  the  statistical  parameters  are  very  useful  for  identifying  the 
nuances  of  the  deflection  basin  changes.  The  linear  changes  in  mean  and  peak  deflection 
reflect  the  manner  in  which  the  vehicle  velocity  influences  the  load  pulses  over  space 
(described  in  Chapter  4).  The  coefficient  of  skewness  shows  that  as  the  speed  increases, 
the  basin  becomes  less  symmetric.  Notice  that  if  the  skewness  curve  were  extended  to 
zero  vehicle  velocity,  it  would  pass  very  nearly  through  zero,  indicating  the  static  basin  is 
perfectly  symmetric. 

7.4.4  Pulse  Duration 

As  described  in  Chapter  4,  the  standard  pulse  duration  of  30  msec  was  selected 
because  it  has  commonly  been  used  by  previous  researchers,  and  because  it  closely 
approximates  the  load  duration  for  an  FWD  test.  Chapter  4  also  showed  that  a  constant 
line  load  could  be  applied  using  any  duration  pulse,  as  long  as  the  amplitude  of  the  pulse 
was  adjusted  accordingly.  This  section  examines  how  the  load  pulse  selection  can 
influence  the  calculated  deflections. 
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Figure  7-27.  Sensitivity  to  Wheel  Speed. 
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Five  different  load  pulses  were  examined,  ranging  between  10  and  50  msec.  Based 
on  Barksdale’s  (1971)  analysis,  a  10  msec  pulse  corresponds  to  a  wheel  velocity  of  over 
75  km/hr  (47  mph),  and  a  wheel  traveling  at  5  km/hr  (3  mph)  would  produce  a  50  msec 
pulse.  However,  the  wheel  velocity  was  held  constant  for  these  analyses  —  only  the  pulse 
duration  was  adjusted. 

Figure  7-28  illustrates  the  computed  deflection  behavior  for  each  pulse  duration. 
The  basin  became  deeper  and  moved  farther  behind  the  wheel  as  the  pulse  duration 
increased. 

The  peak  deflection  consistently  decreased  by  10  pm  for  each  incremental  change 
in  pulse  duration.  The  mean  consistently  increased  (moved  farther  behind  the  wheel)  as 
the  pulse  duration  was  increased  and  the  standard  deviation  increased  the  most  at  the 
longest  pulse  durations.  The  coefficient  of  skewness  increased  almost  linearly  from  0.12 
to  0.47  over  the  range  of  pulse  durations.  The  coefficient  of  kurtosis  was  nearly  constant 
for  the  three  shortest  durations,  and  then  reduced  by  7.5%  as  the  duration  increased  from 
30  to  50  msec,  reflecting  the  flatter  peaks  of  the  longer  duration  basins. 

This  analysis  clearly  shows  that  the  deflection  basin  is  sensitive  to  the  pulse 
duration.  It  seems  that  the  appropriate  duration  should  be  selected  as  a  function  of  the 
vehicle  speed.  However,  when  this  is  done,  the  resulting  deflections  do  not  match  the 
trend  illustrated  in  the  previous  examination  of  wheel  speed.  This  phenomenon  can  be 
attributed  to  the  manner  in  which  the  dynamic  response  due  to  a  load  pulse  is  used  in 
BASIN.  Shorter  duration,  larger  amplitude  pulses  at  the  smallest  radii  have  a  greater 
contribution  to  the  total  deflection  than  longer  duration,  smaller  amplitude  load  pulses 
applied  at  the  same  distances. 
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Figure  7-28.  Sensitivity  to  Load  Pulse  Duration. 
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7.5  Transverse  Deflection  Profile 

All  of  the  previous  results  were  based  on  longitudinal  deflection  basins,  which  show 
the  pavement  deflections  along  the  direction  of  travel.  However,  the  basin  produced  by  a 
moving  wheel  load  is  clearly  three-dimensional.  This  section  examines  the  deflections  in 
the  transverse  direction  (perpendicular  to  the  direction  of  travel)  for  one  parameter,  base 
modulus.  The  deflection  basin  at  a  transverse  plane  intersecting  the  wheel  path  0.2  m 
(7.9  in)  behind  the  wheel  was  calculated  for  each  value  of  the  base  modulus.  This  plane 
was  selected  to  ensure  the  peak  deflection  for  the  entire  basin  was  included. 

The  base  moduli  were  varied  over  the  same  range  used  in  the  longitudinal  analysis, 
from  5  to  50  ksi  (35  to  345  MPa).  Figure  7-29  illustrates  the  differences  in  the 
longitudinal  and  transverse  deflection  profiles. 

The  basins  were  perfectly  symmetric  about  the  peak  deflection,  as  would  be 
expected  for  the  transverse  profile.  The  peak  deflection  and  width  of  the  basins  appeared 
to  be  nearly  identical  to  the  longitudinal  profiles. 

As  illustrated  in  Figure  7-30,  the  peak  deflection  increased  by  20%  when  the  base 
modulus  was  reduced  from  the  standard  to  the  minimum.  As  the  base  modulus  increased 
from  the  standard  to  the  maximum,  the  peak  deflection  decreased  almost  linearly,  of  2 
pm  for  every  1  ksi  increase  in  modulus.  The  mean  was  essentially  zero  for  all  base 
moduli,  and  the  standard  deviation  steadily  increased  as  the  base  modulus  increased.  The 
coefficient  of  skewness  was  nearly  zero  for  all  base  moduli,  confirming  that  the  basins 
were  symmetrical  about  the  mean.  The  coefficient  of  kurtosis  decreased  in  a  non-linear 
manner.  When  the  modulus  was  decreased  from  the  standard  to  the  minimum,  the 
kurtosis  increased  20%,  but  when  the  modulus  was  increased  from  the  standard  by  an 
equal  increment,  the  kurtosis  decreased  by  only  11%. 
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Figure  7-29.  Longitudinal  and  Transverse  Deflection  Basins. 
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Figure  7-30.  Transverse  Profile  --  Sensitivity  to  Base  Modulus. 


The  transverse  basins  behave  similarly  to  the  longitudinal  basins,  except  that  the 
skewness  is  removed.  The  peak  deflections,  standard  deviations,  and  coefficient  of 
kurtosis  are  nearly  identical  for  the  transverse  and  longitudinal  profiles  —  the  small 
differences  are  likely  caused  by  the  skewness  of  the  longitudinal  basins. 


Chapter  8 

SUMMARY,  CONCLUSIONS,  AND  RECOMMENDATIONS 
8.1  Summary 

This  thesis  presented  the  background,  conceptual  development,  and  implementation 
of  a  mathematical  model  for  predicting  pavement  deflections  caused  by  a  moving  wheel 
load.  The  development  of  this  model  was  described  in  the  preceding  chapters. 

Chapter  2  described  a  variety  of  existing  nondestructive  tests  that  are  used  for 
pavement  analysis.  Static,  dynamic,  surface  wave,  and  rolling  wheel  testing  devices  and 
their  historical  applications  were  reviewed. 

Chapter  3  presented  several  types  of  pavement  deflection  models.  First,  the 
simplest  models,  static  linear  elastic,  were  described.  Next,  models  that  can  characterize 
the  pavement  system  as  composed  of  non-linear  elastic  or  viscoelastic  materials  were 
described.  Finally,  dynamic  models  used  for  predicting  transient  stationary  and  moving 
load  deflections  were  explained. 

In  Chapter  4,  the  Green’s  function  method  of  solution  was  applied  to  dynamic 
layered  elastic  analysis.  The  Green’s  functions  used  in  this  model  computed  the  vertical 
deflections  at  the  surface  of  a  pavement  system  resulting  from  a  unit  harmonic  load 
applied  over  a  circular  area.  The  calculated  response  was  a  function  of  the  frequency, 
layer  material  properties,  and  distance  from  the  applied  load.  Fourier  superposition 
analysis  in  combination  with  the  Green’s  function  solution  was  used  to  predict  the 
dynamic  pavement  deflection  response.  Details  of  the  loading  representation  for  a  FWD 
impulse  and  rolling  wheel  were  described. 

The  conceptual  development  of  the  computer  program  BASIN  was  described  in 
Chapter  5.  For  convenience,  a  simulated  deflection  response  function  was  initially  used 
while  the  bookkeeping  methods  were  developed.  This  initial  development  described  the 
modeling  of  five  increasingly  complex  load  and  observation  scenarios.  The  last  portion 
of  the  chapter  described  the  modifications  required  to  insert  the  complex- valued  Green’s 
functions  into  the  program  BASIN. 
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The  results  of  the  deflection  model  validation  were  presented  in  Chapter  6.  The 
rolling  load  model  was  modified  to  determine  dynamic  deflections  caused  by  a  falling 
weight  deflectometer  so  that  predicted  and  measured  deflections  could  be  directly 
compared. 

The  results  of  a  detailed  parametric  study  of  BASIN  were  presented  in  Chapter  7. 
Significant  input  variables  were  systematically  adjusted  to  determine  the  sensitivity  of  the 
deflection  basin  to  each  variable.  Both  layer  material  properties  and  loading  parameters 
were  investigated. 

8.2  Conclusions 

A  number  of  conclusions  were  developed  from  the  results  of  the  parametric  study 
using  the  new  pavement  deflection  model,  coded  in  the  computer  program  BASIN.  The 
conclusions  are  organized  according  to  their  variable  category. 

Algorithm  Considerations 

Algorithm  variables  were  examined  to  determine  the  values  that  produced  the  best 
accuracy  with  the  least  computation  time.  Several  significant  results  are  listed  below. 

1 .  The  thickness  of  the  sublayers  used  to  discretize  the  pavement  must  be 
sufficiently  thin  to  produce  an  accurate  result.  The  parametric  study  showed 
that  the  solution  converges  to  a  constant  deflection  basin  shape,  beyond  which 
thinner  layers  will  not  produce  a  more  accurate  result. 

2.  The  deflection  basin  shape  does  not  change  appreciably  when  the  cutoff 
frequency  used  to  characterize  the  load  is  larger  than  95  Hz. 

3.  The  spacing  between  individual  load  pulses  must  be  less  than  or  equal  to  1 
msec  to  guarantee  an  accurate  solution. 

Material  Considerations 

The  material  variables  were  studied  to  determine  how  their  values  changed  the 
basin  shape.  The  sensitivity  of  the  basin  shape  to  each  variable  may  provide  clues  about 
which  material  parameters  can  be  backcalculated. 
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4.  The  shape  of  the  predicted  deflection  basin  was  insensitive  to  material 
density.  None  of  the  statistical  parameters  were  changed  by  increasing  or 
decreasing  the  density  from  the  standard  case. 

5.  The  deflection  basin  shape  was  relatively  insensitive  to  the  Poisson’s  ratio  of 
each  layer.  Of  the  statistical  parameters,  skewness  and  kurtosis  were  affected 
most  by  the  Poisson’s  ratio,  while  the  peak  deflection  and  standard  deviation 
were  virtually  independent  of  the  Poisson’s  ratio. 

6.  The  deflection  basin  shape  was  very  sensitive  to  the  AC  modulus.  As 
expected,  the  peak  deflection  changed  more  dramatically  than  the  other 
indicators.  Surprisingly,  the  coefficient  of  kurtosis  was  more  sensitive  to  the 
AC  modulus  than  the  standard  deviation.  As  a  result,  there  may  be  more 
potential  for  correlating  AC  modulus  to  kurtosis  than  standard  deviation.  This 
can  be  explained  by  the  manner  in  which  the  basins  changed  —  almost  all  of 
the  change  occurred  very  close  to  the  peak  deflection.  Therefore,  the  spread 
of  the  basin  was  not  changed  as  much  as  its  peak. 

7.  The  deflection  basin  shape  was  very  sensitive  to  the  base  modulus.  Like  the 
AC  modulus,  the  coefficient  of  kurtosis  appeared  to  be  more  sensitive  to  the 
base  modulus  than  the  standard  deviation. 

8.  The  deflection  basin  shape  was  extremely  sensitive  to  the  subgrade  modulus. 
Both  the  peak  and  the  width  of  the  basin  changed  significantly  with  the 
subgrade  modulus.  In  this  case,  the  standard  deviation  was  more  sensitive  to 
the  subgrade  modulus  than  the  other  statistical  parameters. 

9.  The  deflection  basin  shape  was  extremely  sensitive  to  the  stiff  layer  modulus. 
Unlike  the  other  moduli  tested,  the  standard  deviation  of  the  basin  rapidly 
decreased  as  the  stiff  layer  modulus  increased.  The  coefficients  of  skewness 
and  kurtosis  were  also  very  sensitive  to  the  stiff  layer  modulus. 

10.  The  deflection  basin  changed  dramatically  as  the  AC  thickness  changed.  Peak 
deflections  increased  rapidly  as  the  AC  became  thinner.  The  standard 
deviation  decreased  linearly  as  the  AC  thickness  decreased.  Very  thick  AC 
layers  (greater  than  300  mm)  resulted  in  very  broad  deflection  basins  with 
small  dips  near  the  peak  deflection. 

11.  The  deflection  basin  shape  was  very  insensitive  to  the  thickness  of  the  base 
layer.  Almost  no  changes  in  the  statistical  parameters  or  the  general  shape  of 
the  basin  were  observed  for  different  base  layer  thicknesses. 

12.  The  deflection  basin  shape  was  somewhat  sensitive  to  the  subgrade  thickness. 
Unlike  the  behavior  determined  for  the  AC  and  base  thickness,  the  peak 
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deflection  increased  as  the  subgrade  thickness  increased.  This  behavior 
can  be  explained  by  examining  the  stiffness  of  each  layer.  The  subgrade  is 
softer  than  all  the  other  layers.  Therefore,  it  logically  follows  that  a  thicker 
soft  layer  would  result  in  larger  deflections. 

13.  Only  very  small  changes  in  the  deflection  basin  result  from  combined 
variations  of  the  thickness  and  modulus  of  the  AC  or  base  layers.  The 
variations  in  moduli  were  backcalculated  using  static  models,  but 
imexpectedly,  the  dynamic  analysis  used  by  BASIN  did  not  compute 
significantly  different  basins.  These  results  suggest  that  attempts  to 
simultaneously  backcalculate  moduli  and  layer  thickness  will  require  very 
accurate  measurement  of  the  deepest  part  of  the  deflection  basin. 

Loading  Considerations 

Loading  variables  were  studied  to  determine  their  effects  on  the  deflection  basin 
shape.  The  sensitivity  of  the  deflection  basin  shape  to  each  of  these  parameters  may 
provide  clues  about  the  optimum  operating  conditions  for  a  rolling  wheel  deflectometer. 


14.  The  deflection  basin  peak  was  linearly  related  to  the  applied  wheel  load.  The 
other  statistical  parameters  were  constant  for  all  applied  loads.  This  is  a  direct 
result  of  the  assumed  linear-elastic  stress-strain  behavior  of  the  pavement 
materials. 

15.  Only  the  deepest  portion  of  the  deflection  basin  was  sensitive  to  the  tire 
pressure.  Low  tire  pressures  resulted  in  lower  peak  deflections,  but  outside 
the  tire  contact  zone,  the  basin  was  virtually  identical  to  that  produced  by 
higher  pressure  tires. 

16.  The  wheel  speed  has  a  significant  effect  on  the  deflection  basin.  Both  the 
peak  deflection  and  the  mean  were  directly  proportional  to  the  speed.  In  fact, 
wheel  speed  was  the  only  parameter  that  influenced  the  mean.  The  skewness 
also  was  significantly  affected  by  wheel  speed.  At  creep  speeds,  the  basin  was 
nearly  symmetric;  however,  as  the  speed  increased  the  skewness  of  the  basin 
became  positive. 

17.  The  pulse  duration  selected  for  the  analysis  has  a  significant  impact  on  the 
shape  and  location  of  the  deflection  basin.  Short  pulses  produced  deeper, 
more-symmetrical  basins  with  mean  values  closer  to  the  wheel.  Long  pulses 
have  the  opposite  effect.  The  method  in  which  the  relationship  between  load 
pulses  duration,  load  pulse  spacing,  and  wheel  speed  should  be  applied  in 
BASIN  is  not  yet  clearly  understood.  Computations  completed  using  load 
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pulses  corresponding  to  the  wheel  speed  produced  illogical,  inconclusive 
results. 


The  final  analysis  of  the  parametric  study  compared  longitudinal  to  transverse 
deflection  profiles.  These  results  may  begin  to  show  the  advantages  of  one  measurement 
profile  over  the  other. 

18.  Transverse  deflection  profiles  can  show  similar  trends  to  longitudinal  profiles; 
however,  the  skewness  of  the  function  is  removed  in  the  transverse  profile  -- 

i.e.  the  basin  is  perfectly  symmetric.  Transverse  profiles  may  be  better  for 
backcalculation  because  the  absence  of  changes  in  skewness  would  help 
clarify  trends  in  the  other  statistical  parameters.  However,  if  transverse 
profiles  are  measured,  the  distance  between  the  transverse  profile  and  the  peak 
longitudinal  deflection  must  be  accurately  determined  so  that  the  relative 
depth  of  deflection  is  understood.  Also,  parameters  that  produce  changes  in 
skewness  would  not  be  backcalculated  as  accurately  because  the  transverse 
profiles  cannot  recognize  these  changes. 


8.3  Recommendations  for  Further  Study 

The  forward  model  for  predicting  pavement  deflections  caused  by  a  moving  wheel 
load  is  just  one  component  of  a  much  larger  problem.  For  the  model  to  be  of  practical 
value  for  backcalculating  pavement  properties,  several  areas  require  further  study. 

1 .  The  relationship  between  load  pulse  duration,  vehicle  speed,  and  load  pulse 
spacing  should  be  further  investigated.  A  clear  imderstanding  of  these 
parameters  is  required  to  implement  the  model  properly. 

2.  Only  a  limited  number  of  deflection  sensors  can  be  used  to  measure  pavement 
deflections.  Therefore,  it  is  important  to  maximize  the  information  obtained 
fi*om  each  sensor.  The  optimum  location  of  each  sensor  should  be  determined 
based  on  the  sensitivity  of  the  deflection  basin  to  expected  material  parameters 
and  the  most  probable  operating  speed  and  load.  If  possible,  sensors  should 
be  located  to  allow  an  accurate  definition  of  the  peak  deflection  zone,  because 
this  area  may  be  most  sensitive  to  thickness  and  modulus. 

3.  The  statistical  parameters  used  to  characterize  the  deflection  basins  in  the 
parametric  study  may  prove  useful  for  very  rapid  direct  correlations  with  layer 
modulus  and  thickness.  The  application  of  these  parameters  to  FWD 
deflection  basins  may  help  develop  new  correlations  for  stationary  tests  too. 
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4.  The  Green’s  function  solution  developed  by  Kausel  and  Peek  (1982)  is 
capable  of  calculating  vertical  or  horizontal  displacements  within  the 
pavement  structure.  The  program  BASIN  could  be  modified  to  calculate 
deflections  at  critical  locations  in  the  pavement  structure  (i.e.  horizontal 
deflections  at  the  bottom  of  the  AC  layer).  From  this  information,  dynamic 
pavement  damage  models  could  be  developed  to  predict  the  life  expectancy  of 
a  pavement  structure. 

5.  The  FORTRAN  program  BASIN  required  about  70  minutes  to  compute 
deflections  at  30  locations  using  a  486/66  MHz  microcomputer  with  8  Mb 
RAM.  This  computation  time  is  much  too  long  for  practical  application  of  the 
model.  However,  the  modified  program  FWDDROP  can  calculate  the  time 
history  deflection  for  seven  sensor  locations  in  only  1.5  minutes.  Therefore, 
this  model  may  have  some  current  practical  value  for  backcalculation.  As 
computer  processing  speeds  increase,  the  required  computation  times  'will 
undoubtedly  decrease. 

6.  When  deflection  measurements  from  rolling  weight  deflectometers  are 
available,  the  accuracy  of  the  model  should  be  verified. 
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Appendix  A 
BASIN  CODE 


A.1  Primary  BASIN  Modules 

BASIN  is  comprised  of  a  primary  module  and  tliree  main  subroutines.  The 
primary  module  controls  the  order  of  the  program  execution  and  performs  basic 
“bookkeeping”  functions.  The  LOADDAT  subroutine  reads  the  loading  and  velocity 
information  recorded  by  the  rolling  weight  deflectometer.  Also,  avetime,  timeback,  and 
other  parameters  used  to  characterize  individual  load  pulses  are  loaded  into  the  program. 
The  subroutine  GREEN  was  generously  provided  by  Professor  Roger  Meier.  This 
subroutine  calculates  the  Green’s  function  response  for  all  frequencies  within  a  specified 
bandwidth.  LOADTIME  is  called  by  GREEN  after  the  response  at  eaeh  frequency  is 
known.  Subroutine  LOADTIME  performs  Fourier  superposition  analysis  to  create  the 
time  history  of  deflection  for  each  load  pulse.  Within  LOADTIME,  the  subroutine  FFT 
is  called  to  complete  forward  and  inverse  Fast  Fourier  Transforms.  The  subroutine  FFT 
was  taken  directly  from  Press  et  al.  (1989).  Finally,  LOADTIME  selects  the  proper 
deflection  from  each  time  history  and  calculates  the  total  deflection  at  each  deflection 
sensor  location. 

A.2  BASIN  Input  and  Output 

BASIN  uses  four  input  files. 

pulse.dat  This  file  contains  the  time  increment,  avetime,  timeback, 
pulse  duration,  total  load  amplitude,  and  wheel  speed. 

rwd.dat  This  file  contains  a  simulated  list  of  output  data  measured 
by  a  rolling  weight  deflectometer.  rwd.dat  is  created  using 
the  program  CRE8DAT. 
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input  dat  This  file  is  used  by  GREEN  and  contains  the  material 

layer  properties,  loaded  radius,  and  frequency  information. 
See  the  subroutine  data  in  GREEN  for  formatting  details. 

sense.dat  The  initial  location  of  the  wheel  and  the  locations  of  all 
deflection  sensors  are  listed  in  this  file. 


BASIN  writes  the  output  to  one  file  -  output  out  The  output  file  first  lists  the 
basin  number  and  the  time  and  location  of  the  wheel  when  the  calculation  was  made. 
Next,  the  coordinates  for  each  deflection  sensor  and  the  deflection  at  that  sensor  is  listed. 


A.3  BASIN  Code 


program  BASIN 
c 

c  - 

c  This  program  computes  deflection  resulting  from  a  rolling  wheel  load, 
c  BASESf  uses  Green’s  functions  that  determine  the  steady-state  response 
c  for  a  unit  vibrating  disk  load.  The  green's  functions  are  computed 
c  in  the  subroutine  GREEN,  written  by  Roger  Meier,  with  minor  modifications 
c  for  this  purpose. 

c  NOTE:  When  compiling  on  MS  FORTRAN  Powerstation,  two  warnings  will  be 
c  given  that  could  not  be  removed.  However,  they  do  not  prevent 
c  program  execution.  They  are  both  in  line  1059  and  read: 
c  warning  F4016:  LOVE  :  formal  argument  E  :  type  mismatch 
c  warning  F4016:  LOVE  :  formal  argument  V  :  type  mismatch 

c  - ^ - 

implicit  real*  8  (a-h,0“Z) 
parameter  (iout=16,isense=25) 

parameter  (ntimemax=20 1  ,nmbkmax=20 1  ,nbasnmax=ntimemax-nmbkmax+ 1 , 

1  nterms= 1 024,ndefinax=40) 

parameter(itim=86) 
integer*2  ihr,imin,isec,ihun 

real*  8  dt,time(ntimemax), velo(ntimemax),pleft(ntimemax), 

1  pright(ntimemax) 

real*  8  def  1  (ntimemax),def2(ntimemax),def3(ntimemax) 

real*  8  xinit,yinit 

real*  8  xd(ndefrnax),yd(ndefmax) 

real*  8  xwheel(ntimemax) 

real*  8  xobserv(nbasnmax,ndefinax) 

real*  8  yobserv(ndefmax) 

real*  8  radius(nmbkmax,ndefinax) 
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real*  8  findef(nbasnmax,ndefmax),calcdef(ndefmax) 
real*  8  avetime,  pulseamp,totalamp 

c - Start  timer 

open  (itim,file-timer.out', status- new’) 
call  gettim(ihr,imin,isec,ihun) 
write(itim,*)  ihr,':',imin,':',isec 

open  (isense,file-sense,dat', status- old') 
open  (ioutjfile  =  'output.ouf,  status  =  'new') 

c - Load  input  data  from  RWD 

call  LOADAT  (dt,ntime,avetime,timeback,numback,nbasin, 

1  fwda,fwdb,totalamp,pulseamp, 

2  time,velo,pleft,pright, 

3  defl,def2,def3, 

4  ntimemax) 
c 

c - Read  in  basic  location  paramaters 

c  xinit  =  initial  x  position  of  wheel 
c  yinit  =  initial  y  position  of  wheel 
c  xd(i)  =  x  position  of  deflection  sensor  i 

c  yd(i)  =  y  position  of  deflection  sensor  i 

c - 


read(isense,*)  xinit, yinit,ndef 
do  i=l,ndef 

read(isense,*)  xd(i),yd(i) 
end  do 
close(isense) 

c - Calculate  the  position  of  the  wheel  at  all  times 

c  Assume  the  wheel  moves  in  the  +x  direction 
c 

xwheel(l)=xinit 
do  i-2,ntime+l 

xwheel(i)=xwheel(i-l)  +  velo(i-l)*avetime 
end  do 

c - Numback  specifies  how  far  backwards  in  time  the  program  should  look  for 

c  calculating  a  response. 

c  Units  for  numback  are  in  terms  of  the  number  of  data  points, 
c 

c - By  default,  the  first  point  at  which  we  can  calculate  a  basin  is  at 

c  xwheel(numback), 
c 

c - Establish  the  observation  point  coordinates  for  all  nbasins. 

do  jdef=l,ndef 
y  observCi  def)=y  d(j  def) 
end  do 

do  i=l,nbasin 

mtemp=numback  +  (i-1) 
do  jdef=l,ndef 

xobserv(i,jdef)=xwheel(mtemp)+xd(jdef) 
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end  do 
end  do 

write(*,*)  ’Observation  point  coordinates  established' 
do  jbasin=l5nbasin 


c - Calculate  matrix  containing  all  radii  [numback  x  ndef] 

do  jdef=l,ndef 

do  iback=l, numback 

radius(iback,jdef)=dsqrt((xwheel(iback+jbasin- 1  )- 

1  xobserv(jbasinJdef))**2  + 

2  yobserv(jdef)**2) 
end  do 

end  do 


write(*,*)  'Radius  matrix  established' 


c 

c - Call  GREEN  to  calculate  the  time  history  response 

c  for  each  radii, 
c 

call  GREEN(radius,calcdef,dt, 

1  avetime, numback, ndef, 

2  fwda,fwdb,pulseamp) 

write(*,*)  'THRs  complete' 

c - Transfer  the  calcdef  variable  to  a  permanent  fmdef  variable 

do  jdef=l,ndef 

findef(jbasinjdef)=calcdef(jdef) 
end  do 

c - Print  the  results  in  the  output  file 


Avrite(iout,*)  'Basin  #'Jbasm 

write(iout,'(2e  1 4.6)')  time(jbasin- 1  +numback+ 1 ), 

1  xwheel(jbasin- 1  +numback+ 1 ) 

do  jdef=l,ndef 

write(iout,’(3el4.6)')  xd(jdef),yd(jdef),fmdef(jbasin,jdef) 
end  do 

write(iout,*) '  ' 

write(*,*)  'Basin  #'Jbasin,'  complete' 
end  do 

c - Stop  timer 

call  gettim(ihr,imin,isec,ihun) 
write(itim,*)  ihr,':',imin,':',isec 
close(itim) 
close(iout) 

stop 

end 


c 


subroutine  LOAD  AT  (dt,ntime, avetime,timeback, numback, nbasin, 
1  fwda,fwdb,totalamp,pulseamp. 


2 

3 

4 


time,velo,pleft,pright, 

defl,def2,def3, 

ntimemax) 


c - 

c  Subroutine  LOAD  AT 

c  This  routine  first  loads  several  basic  parameters  used  throughout 
c  the  main  program.  These  are  in  2  input  data  files/ 
c 

c  pulse.dat 
c  dt  =  time  increment 

c  avetime  =  time  spacing  between  load  pulses 

c  timeback  =  span  of  time  used  to  determine  the  load  pulses  used 

c  fwda  =  start  of  load  pulse  (should  always  be  zero) 

c  fwdb  =  end  of  load  pulse 

c  totalamp  =  total  line  load  applied  by  the  wheel 

c 

c  Next,  This  routine  loads  data  from  the  RWD  into  the  program, 
c  For  model  development,  this  data  wascreated  by  using  the 
c  program  CRE8DAT. 

c  Data  should  be  in  seven  columns  of  real  numbers  in  the  order: 
c 

c  rwd.dat 

c  Time  Velocity  Left  Wheel  Right  Wheel  Deflect  Deflect  Deflect 
c  Load  Load  1  2  3 

c 

c  Total  number  of  terms  to  be  loaded  is  specified  by  ntime. 

c - 

c 

implicit  real*  8  (a-h,o-z) 
parameter(idat=  1 4,loadinp=  1 0) 

real*  8  dt,time(ntimemax), velo(ntimemax),pleft(ntimemax), 

1  pright(ntimemax) 

real*  8  def  1  (ntimemax), def2(ntimemax),def3  (ntimemax) 
real*  8  avetime, totalamp,pulseamp 
real*8  dummy 

open(loadinp,file=‘pulse.dat’,status='old') 
read(loadinp,*)  dt, avetime, timeback 
read(loadinp,*)  fwda, fwdb, totalamp 

ntime=timeback/avetime 
numback  =  timeback/avetime 
dummy=(fwdb-fwda)/2 
pulseamp  =  totalamp/(dummy/avetime) 
nbasin  =  ntime-numback+1 
c 

open(idat, file- rwd.dat’, status=’old’) 
do  i=l,ntime+l 

read  (idat,*)  time(i),velo(i),pleft(i),pright(i), 

1  defl(i),  def2(i),  def3(i) 

end  do 
c 


close(loadinp) 


close(idat) 

return 

end 


c 

c 

c 

c 

c 

subroutine  LOADTIME(dt,uz,nfreq,calcdef, 

1  avetime,numback,ndef, 

2  fwda,fwdb,pulseamp) 


c  This  routine  creates  an  impulse  load  shape  and  converts  it  to 
c  the  frequency  domain  using  an  FFT.  Next,  it  multiplies  the 

c  load  times  the  green’s  function  response  for  each  frequency, 

c  Finally,  an  inverse  FFT  of  the  result  is  used  to  produce  a 
c  time  history  response  (thr). 


parameter(nterms= 1 024,nterms2=2  *nterms) 
parameter(nmbkmax=20 1  ,ndefrnax=40,maxfr=40) 
parameter(pi=3. 14 1 5926535898) 
c  parameter(ibug=86) 
implicit  real* 8  (a-h,o-z) 
real*  8  fwda,fwdb,pulseamp 
real*  8  dt,fwdtime 
real*  8  avetime 

real*  8  pload(nterms),pdata(nterms2),data(nterms2) 
real*  8  fdata(nterms2),tdata(nterms2),rdata(nterms2) 
real*  8  thr(nterms) 

real*  8  tdef(ndefmax),deflec(nmbkmax),calcdef(ndefmax) 
complex*  16  uz(maxfr,nmbkmax,ndefhiax) 
complex*  16  fload(nterms),truncp(nterms),eye 
complex*  16  fresp(nterms),gmresp(nterms),static 

eye=cdsqrt(-l) 

c - Create  loading  shape  in  time  domain  based  on  input  parameters 

do  i=l,nterms 
fwdtime=(i”l)*dt 

if(fwdtime.gt.fwda.and.fwdtime.lt.fwdb)  then 
pload(i)=0.5*pulseamp*cos((fwdtime-fwda)*2*pi/ 

1  (fwdb-fwda)+pi)+pulseamp*0.5 

else 

pload(i)=0.0 
endif 
end  do 

c - Stack  the  string  of  numbers  into  a  single  string  of  alternating 

c  real  and  imaginary  values 
do  i=l,nterms 
pdata(  1 +(i- 1  )*  2)=pload(i) 


pdata(2+(i- 1)*2)=0.0 
end  do 

c - Convert  to  frequency  domain  using  FFT  function 

isign=l 

data=pdata 

call  FFT(data,nterms,isign) 
fdata=data 

c/RGT  Added  to  fix  reversed  sign  of  phase  angle 
do  i=l,nterms 

fdata(2+(i-l)*2)=-fdata(2+(i-l)*2) 
end  do 
c/RGT 

c - Reassemble  fload  into  complex  number  format 

do  i=l,nterms 

fload(i)=fdata(  1  +(i- 1  )*2)+eye*  fdata(2+(i"  1  )*2) 
end  do 

c - Truncate  the  frequencies  to  use  fewer  to  define  the  load 

write(*,*)  'nfreq=’,nfreq 
do  i=l5nfreq 
truncp(i)=fload(i) 
end  do 

do  i=nfreq+ 1  ,nterms-(nfreq- 1 ) 
truncp(i)=0 
end  do 

do  i=nterms,nterms-(nfreq-2),-l 
truncp(i)=fload(i) 
end  do 

write(*,*)  ’Entering  main  LOADTIME  loop’ 
c - Convert  uz  into  double-sided,  green's  function  string 

do  jdef=l,ndef 
tdef(jdef)=0.0d0 
do  iback=  1  ,numback 
static=uz(l  ,iback  jdef) 
gmresp(  1  )=dreal(static) 
do  i=2,nfreq 

gmresp(i)=uz(i,ibackjdef) 
end  do 

do  i=nfreq+ 1  ,(nterms-(nfreq- 1 )) 
gmresp(i)=0.0d0 
end  do 
ncount=2 

do  i=nterms,(nterms-(nfreq-2)),-l 
ntemp=i-nterms+ncount 
gmresp(i)=dreal(uz(ntemp,iback,jdef))- 
1  eye*dimag(uz(ntemp,iback,jdef)) 

ncount=ncount+2 
end  do 


c - Multiply  the  green's  function  by  the  load  in  the  frequency  domain 


do  i=l,ntenns 

fresp(i)=truncp(i)  *  gmresp(i) 
end  do 

c - Restack  the  response  data  into  proper  format  for  inverse  FFT  evaluation 

do  i=l,nterms 

rdata(  1  +(i- 1  )*2)=dreal(ff  esp(i)) 
rdata(2+(i- 1  )*2)=-dimag(fresp(i)) 
end  do 

c - Convert  back  to  the  time  domain  using  an  inverse  FFT 

isign=-l 

data=rdata 

call  FFT(data,nterms,isign) 
tdata=data 

c - Unstack  tdata  into  a  string  of  real  time  values,  thr, 

c  neglecting  the  imaginary  part  of  the  string, 
do  i=l,nterms 
thr(i)=tdata(  1  +(i- 1 )  *  2) 
end  do 

c/RGT— If  thr  values  at  small  and  large  times  are  less  than  zero,  set  to  zero. 
ntrip=0 
do  i=l,512 
if(ntrip.eq.O)  then 
if(thr(i).lt.0.0)  then 
thr(i)-0.0d0 
else 
ntrip=l 
end  if 
else 

continue 
end  if 
end  do 

do  i=513,nterms 
thr(i)=0.0d0 
end  do 
c/RGT - 


c -  Select  the  proper  deflections  and  sum  them 

ntemp=(numback-iback+ 1  )*  (avetime/dt)+ 1 

deflec(iback)=thr(ntemp) 

tdef(jdef)=tdef(jdef)+deflec(iback) 

c - close  the  numback  loop 

end  do 

calcdef(jdef)=tdef(jdef) 

c - close  the  ndef  loop 

end  do 

c  close(ibug) 
return 


end 
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subroutine  FFT(data,nterms,isign) 

c - 

c  This  function  takes  the  FFT  or  IFFT  of  the  string  of  complex 
c  numbers,  data,  which  contains  nn  complex  terms.  If  isign  1, 
c  and  FFT  is  performed,  if  isign  =  -1,  an  IFFT  is  performed, 
c  Data  must  be  packed  such  that  data(l)  is  the  real  portion 
c  of  the  first  input  term  and  data(2)  is  the  imaginary  portion, 
c  Copied  from  "Numerical  Recipes"  by  ??? 


integer*4  isign,nterms 

real*  8  data(nterms*2) 

integer*4  i,istepj,m,mmax,n 

real*  8  tempi, tempr 

real*  8  theta, wi,wpi,wpr,wr,wtemp 

n=2*nterms 

j=l 

do  i=l,n,2 
if(j.gt.i)then 
tempr=data(j) 
tempi=data(j+l) 
data(j)=data(i) 
data(j+ 1  )=data(i+ 1 ) 
data(i)=tempr 
data(i+l)=tempi 
end  if 
m=n/2 

1  if((m.ge.2).and.(j.gt.m))  then 

m=m/2 
goto  1 
end  if 
j=j+m 
end  do 
mmax=2 

2  if  (n.gt.mmax)  then 

istep=2*mmax 

theta=6.2831 85307 17959d0/(isign*mmax) 
wpr=-2.d0*sin(0.5d0*theta)**2 
wpi=sin(theta) 
wr=l.d0 
wi=0.d0 
do  m=l,mmax,2 
do  i=m,n,istep 
j=i+mmax 

tempr=sngl(wr)*data(j)-sngl(wi)*data(j+ 1 ) 
tempi=sngl(wr)*data(j+l)+sngl(wi)*data(j) 
dataG)=data(i)-tempr 
data(j+ 1  )=data(i+ 1  )-tempi 


data(i)=clata(i)+tempr 
data(i+ 1  )=data(i+ 1  )+tempi 
end  do 
wtemp=wr 

wr=wr*wpr-wi*wpi+wr 
wi=(wi*wpr)+(wtemp*wpi)+wi 
end  do 
mmax=istep 
goto  2 
end  if 

c/RGT  Add  normalizing  term  for  inverse  transform 
if(isign.eq.-l)  then 
do  i=l,nterms 

data(  l+(i- 1  )*2)=data(  1  +(i"  1  )*2)/nterms 
end  do 
else 

continue 
end  if 
c/RGT 

return 

end 


subroutine  GREEN(radius,calcdef,dt, 

1  avetime,numback,ndef, 

2  fwda,fwdb,pulseamp) 
c 

c  - 

c/RGT 

c  Obtained  from  Roger  Meier 
c/RGT 

c  this  is  a  modified  version  of  ignacio  sanchez-salinaro's  program 
c  to  find  the  green's  function  solutions  to  elastodynamic  problems 

c  - 

c  -  the  program  and  subroutine  banners  have  been  cleaned  up 

c  -  the  file  code  it=17  has  been  moved  to  the  parameter  statement 

c  -  the  temporary  files  have  been  eliminated  (they're  never  used) 
c  -  the  inputs  for  the  load  have  been  moved  to  subroutine  data 
c  -  the  fixed  file  code  it  replaces  the  variable  file  code  loadt 

c  -  fi'equency  has  been  added  to  the  output  to  file  code  io 

c  -  frequency  has  been  removed  from  the  output  to  file  code  it 

c  -  all  output  moved  from  the  subroutines  into  the  main  program 
c  -  added  close  statements  for  all  of  the  input/output  file  codes 
c  -  subroutine  hankel  replaced  with  eduardo  kausel's  version 

c  -  subroutine  jbessel  replaced  with  eduardo  kausel's  version 

c  -  added  code  to  automatically  create  the  nodes  w/in  each  layer 

c  -  added  array  ns  to  specify  the  number  of  segments  per  layer 

c  -  added  an  exhaustive  input  data  definition  section 


-  collapsed  the  input  read  statements  into  implicit  do-loops 

-  removed  references  to  loading  cases  that  are  not  supported 

-  added  loads  to  the  displacement  calculations  for  line  loads 

-  reformatted  subroutine  displac  to  enhance  its  readability 

-  moved  all  i/o  from  remaining  subroutines  to  main  program 

-  removed  all  debugging  i/o  that  had  been  commented  out 

-  added  comment  statements  to  subroutines  data  and  echo 

-  removed  subroutine  line  (which  was  never  called  anywhere) 

-  moved  four  and  rr  into  subroutine  boupra  as  parameters 

-  eliminated  cs  and  cp  arrays  (never  used  except  for  halfspace) 

-  moved  cs  and  cp  calculations  for  halfspace  into  hspace 

-  moved  call  to  subroutine  boupra  inside  frequency  do-loop 

-  eliminated  temporary  variables  used  inside  frequency  do-loop 


see  subroutine  data  for  a  description  of  the  required  input  data 


implicit  real*  8  (a-h,o-z) 

parameter  (nterms=1024,mlay=5,msub=40,statfr=0.01, 

1  mdim=msub+ 1  ,mdim2=2*mdim,mdim3=3  *mdim, 

2  pi=3. 1415926535898, pi2=6.2831853071796) 
parameter(nmbkmax=20 1  ,ndefmax=40,maxfr=40) 
parameter(in=12) 

real*  8  rock 

integer*4  layer(mdim),jfreq 
mteger*4  iwk(mdim3) 
real*  8  h(mdim) 

real*  8  th(mlay),ns(mlay),  vs(mlay),rro(mlay),anu(mlay),damp(mlay) 
real*  8  mil  (mdim),ml2(mdim),mr  1  (mdim2),mr2(mdim2) 
real*  8  load 

real*  8  dt,radius(nmbkmax,ndefmax) 
real*  8  avetime,fwda,fwdb,pulseamp 
real*  8  calcdef(ndefmax) 
real*  8  holdofr 
complex*  16  lame(mlay,3) 

complex*  16  all(mdim),al2(mdim),gll(mdim),gl2(mdim) 
complex*  16  dll(mdim),drl(mdim),dr2(mdim),dr3(mdim) 
complex*  16  erl(mdim),er2(mdim),trl(mdim),tr2(mdim) 
complex*  16  qrl(mdim),qr2(mdim),prl(mdim),pr2(mdim) 
complex*  16  arl(mdim2),ar2(mdim2),brl(mdim2),br2(mdim2) 
complex*  16  grl(mdim2),gr2(mdim2) 
complex*  1 6  bl  1  (mdim),bl2(  1 00) 

complex*  16  wkl(mdim3),wk2(mdim2),wk3(mdim2),wk4(mdim2) 
complex*  16  wk5(mdim2),wk6(mdim2),wk7(mdim2),wk8(mdim2) 
complex*  16  crl(mdim2),cr2(mdim2) 
complex*  16  er(mdim2),vr(mdim2,mdim2),wk(mdim2,4) 
complex*  16  uzl,uz(maxfr,nmbkmax,ndefinax) 

open(in,file-input.daf,status=’old') 


call  data  (th,ns,vs,rro,anu,damp,lame, 
1  nfreq,ofr,dfr,ifr,rock, 
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2  r,load,loadt,il,id, 

3  ndist, odist, ddistjidist, 

4  niay, 

5  mlay,in) 

c/RGT-Statfr  is  the  smallest  frequency  calculated  -  used  to  approximate 
c  the  solution  for  a  static  load. 
holdofr=ofr 
ofr=statff 

c/RGT- - 

call  sublay  (layer, h, 

1  th,ns,vs, 

2  nfreq,ofr,dff,ijfr,rock, 

3  nlay,nsub, nodes, nodes2, 

4  mlay,mdim) 


fr  =  off 

dojffeq=  l,nffeq 
om  =  fr*pi2 
oms  =  om*om 

call  boupra  (all,al2,gll,gl2,mll,ml2,arl,ar2,brl,br2, 

1  grl  ,gr2,mrl  ,mr2,dll  ,drl  ,dr2,dr3,erl  ,er2, 

2  tr  1  ,tr2,qr  1  ,qr2,pr  1  ,pr2, 

3  layer, h, 

4  rro,lame, 

5  nsub, 

6  mlay,mdim,mdim2) 

call  hspace  (all,al2,gll,gl2,ml2,arl,ar2,brl,br2, 

1  grl,gr2,mr2,dr2,dr3,er2, 

2  tr2,qr2,pr2, 

3  om,rock, 

4  layer, 

5  lame,vs,anu,damp, 

6  nodes,nodes2, 

7  mlay,mdim,mdim2) 

call  prelove  (bll,bl2,gll,gl2,mll,ml2, 

1  oms, 

2  nodes, 

3  mdim) 

call  love  (bl  1  ,bl2,al  1  ,al2, 

1  el,vl, 

2  wkl  ,wk2,wk3,wk4,wk5, 

3  iwk, nodes, 

4  mdim,mdim2,mdim3 ,99) 
call  prerayl  (crl,cr2,grl,gr2,mrl,mr2, 

1  oms, 

2  nodes, 

3  mdim2) 

call  raylgh  (arl,ar2,brl,br2,crl,cr2, 

1  er,vr, 

2  wk,wkl  ,wk2,wk3,wk4,wk5,wk6,wk7,wk8, 

3  iwk,nodes2, 

4  mdim2,mdim3,99) 
do  jdef=l,ndef 
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c/RGT - Start  loop  to  calculate  response  for  all  radii  at  this  frequency 

do  iback=l,numback 
dist=radius(ibackjdef) 

c/RGT - 

call  displac  (uzl, 

1  er,vr, 

2  r,load,loadt,iUd, 

3  dist, 

4  nodes2, 

5  mdim2,pi) 
uz(jfreq,iback,jdef)=uz  1 

end  do 
end  do 

write(*,*)  'frequency  =',fr, 'complete' 
c/RGT — Revert  back  to  the  original  starting  frequency  when  statfr 
c  is  completed 

if(jfreq.eq,l)  ofr=holdofr 

c/RGT - 

if  (ifr.eq.O)  then 
fr  =  ofr  +  dfr*jfreq 
else 

fr  =  ofr  +  dfr*2**(jfreq-l) 
end  if 

end  do 
c 

c/RGT-Subroutine  to  apply  load  shape  and  do  FFT  conversions 
write(*,*)  'All  frequencies  complete.  Entering  Loadtime' 
call  LOADTIME(dt,uz5nfreq,calcdef, 

1  avetimCjnumbackjndef, 

2  fwda,fwdb,pulseamp) 

C/RGT - 

go  to  100 

99  write  (*,*)  'problems  with  eigenvalues  at  fr  =  ',fr 

100  continue 
close(in) 
return 

end 

c 

c 

c 

c 

c 

subroutine  data  (th,ns,vs,rro,anu,damp,lame, 

1  nfreq, ofr, dfr,ifr, rock, 

2  r,load,loadt,il,id, 

3  ndist,odist,ddist,idist, 

4  nlay, 


5 


mlay,in) 
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c  - 

c  input  data 

c 

c  read  (in,*)  nlay,rock 

c  nlay . number  of  material  layers 

c  rock . bottom  boundary  condition 

c  ’r  =  half-space 

c  ’O'  =  rigid 

c  read  (in,*)  (th(i),ns(i),vs(i),rro(i),anu(i),damp(i),i=l,nlay) 

c  th(i) . total  thickness  of  layer  i 

c  ns(i) . number  of  segments  in  layer  i 

c  vs(i) . shear  wave  velocity  in  layer  i 

c  rro(i) . mass  density  in  layer  i 

c  anu(i) . poisson's  ratio  in  layer  i 

c  damp(i) . fraction  of  damping  in  layer  i 

c  read  (in,*)  r,load,loadt,il,id 

c  r . radius  of  applied  load 

c  load . magnitude  of  applied  load 

c  =  force  for  point  loads 

c  =  force  per  length  for  line  loads 

c  =  force  per  area  for  area  loads 

c  loadt . type  of  applied  load: 

c  1  =  point  load  vertical 

c  3  =  disk  load  vertical 

c  il . node  where  load  applied 

c  id . node  where  displacements  output 

c  read  (in,*)  ndist,odist,ddist,idist 

c  ndist . number  of  output  locations 

c  odist . starting  output  location 

c  ddist . output  location  intervals 

c  idist . type  of  progression 

c  0  =  arithmetic 

c  1  =  geometric 

c  read  (in,*)  nfreq,ofr,dfr,ifr 

c  nfreq . number  of  frequencies 

c  ofr . starting  frequency 

c  dfr . frequency  interval 

c  ifr . type  of  progression 

c  0  =  arithmetic 

c  1  =  geometric 

c  - 

implicit  real* 8  (a-h,o-z) 

real*  8  rock 

real*  8  th(mlay),ns(mlay),vs(mlay),rro(mlay),anu(mlay),damp(mlay) 

real*  8  load 

complex*  16  lame(mlay,3) 

complex*  16  gg,ee 

c  read  material  layer  data 


read  (in,*)  nlay,  rock 
write  (*,*)  'read  nlay 


1 
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do  i  =  l^nlay 

read  (in,*)  th(i),ns(i),vs(i),rro(i),anu(i),damp(i) 
gg  =  rro(i)*vs(i)*vs(i)*dcmplx(1.0d0,2.0d0*damp(i)) 
ee  =  2.0d0*gg*(1.0d0+anu(i)) 
lame(i,l)  =  ee 
lame(i,2)  =  gg 

lame(i,3)  =  2.0d0*anu(i)*gg/(1.0d0-2.0d0*anu(i)) 
end  do 

c  read  applied  load  data 

read  (in,*)  r,load,loadt,il,id 

c  read  output  location  data 

read  (in,*)  ndist,odist,ddist,idist 
c  read  frequency  data 

read  (in,*)  nfreq,ofr,dfr,ifr 

return 

end 

subroutine  sublay  (layer, h, 

1  th,ns,vs, 

2  nfreq,ofr,dfr,ifr,rock, 

3  nlay,nsub, nodes, nodes2, 

4  mlay,mdini) 

implicit  real*  8  (a-h,o-z) 
real*  8  rock 
integer*4  layer(mdim) 
real*8  h(mdim) 

real*  8  th(mlay),ns(mlay),vs(mlay) 
c  compute  minimum  and  maximum  frequencies 


fhnin  =  ofr 
if  (ifr.eq.O)  then 
fimax  =  ofr  +  dfr*(nfreq-l) 
else 

frmax  =  ofr  +  dfr*2**(nfreq-2) 
end  if 

print  *, 'minimum  frequency  specified  = ', fhnin 
print  *, 'maximum  frequency  specified  =  ',frmax 

c  compute  minimum  and  maximum  wavelengths 

wavmin  =  wavlen(th,vs,nlay,mlay,frmax) 
wavmax  =  wavlen(th,vs,nlay,mlay,frmin) 

print  *,'minimum  Rayleigh  wavelength  =  ',wavmin 


print  *, 'maximum  Rayleigh  wavelength  =  ’,wavmax 
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c  generate  sublayers  within  each  material  layer 

nsub  =  0 
do  lay  =  l,nlay 
dz  =  th(lay)/ns(lay) 
do  j  =  l,nsOay) 
nsub  =  nsub  +  1 
h(nsub)  =  dz 
layer(nsub)  =  lay 
end  do 
end  do 

c  add  another  sublayer  for  the  transmitting  boundary 

nodes  =  nsub 
if(rock.eq.l)  then 
nodes  =  nodes  +  1 
layer(nodes)  =  nlay 
endif 

nodes2  =  nodes*2 

return 

end 

function  wavlen  (th,vs,nlay,mlay,fr) 
c 

c  computes  approximate  wavelength  of  dispersed  Rayleigh  wave 

c  - 

c 

implicit  real*  8  (a-h,o-z) 
real*  8  th(mlay),vs(mlay) 

w(z)  =  1.0  -  1.36765*exp(-2.2281*z)  +  0.36765*exp(-5.6137*z) 

guess  =  vs(nlay)/fr 
1  above  =  0.0 
depth  =  0.0 
speed  =  0.0 
do  i  =  l,nlay 
depth  =  depth  +  th(i) 
below  =  w(depth/guess) 
speed  =  speed  +  vs(i)*(below-above) 
above  =  below 
end  do 

if  (above.lt.  1.0)  speed  =  speed  +  vs(nlay)*(1.0-above) 

wavlen  ^  speed/fr 

if  (abs(wavlen-guess).lt.0.15)  return 

guess  =  wavlen 

goto  1 


end 
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subroutine 
1 
2 

3 

4 

5 
c 

c  - 

c  computes  displacements  using  Green’s  function 

c  - 

c 

implicit  real*  8  (a-h,o-z) 
real*  8  load 

complex*  16  er(mdim2),vr(mdim2,mdim2) 
complex*  1 6  uzl 
complex*  16  ev,zr,zd 

complex*  1 6  hOd,h  1  d,hOr,h  1  r  j  Od  J 1  d  jOrj  1  r,i  1  r 
complex*  16  vzd,vzl 
complex*  16  pihi 

c  vertical  point  load 

if(loadt.eq.l)  then 
uzl=(0.0d0,0.0d0) 
do  l=l,nodes2 
ev=er(l) 
zd=ev*dist 

call  hankel  (zd,h0d,hld,2) 
vzd=vr(2*id,l) 
vzl=vr(2*il,l) 
uzl  =  uzl  +  vzd*vzl*hOd 
end  do 

uzl  =uzl*load/(0.0d0,4.0d0) 
endif 


c  vertical  disk  load 

if(loadt.eq.3)  then 
uzl=(0.0d0,0.0d0) 
pihi=pi/(0.0d0,2.0d0) 
do  l=l,nodes2 
ev=er(l) 
zr=ev*r 
zd=ev*dist 
if(dist.lt.r)  then 
call  hankel  (zr,h0r,hlr,2) 
call  bessel  (zdJOdjld) 
ilr=pihi*jOd*hlr/ev  -  l.OdO/r/ev/ev 
else 

call  hankel  (zd,h0d,hld,2) 


displac  (uzl, 
er,vr, 

r,load,loadt,il,id, 

dist, 

nodes2, 

mdim2,pi) 


call  bessel  (zrjOrJlr) 

ilr=pihi*jlr*hOd/ev 

endif 

vzd=vr(2*id,l) 
vzl=vr(2*il,l) 
uzl  =  uzl  +  vzd*vzl*ilr 
end  do 

uzl  =uzl*load*r 
endif 

return 

end 


subroutine  boupra  (all,al2,gll,gl2,mll,ml2,arl,ar2,brl,br2, 

1  gr  1  ,gr2,mr  1  ,mr2,dl  1  ,dr  1  ,dr2,dr3  ,er  1  ,er2, 

2  trl,tr2,qrl,qr2,prl,pr2, 

3  layer, 

4  h,rro,lame, 

5  nsub, 

6  mlay,mdim,mdim2) 


generates  boundary  layer  matrices  (a,b,g,m,e,d,n,l,q) 


implicit  real*  8  (a-h,o-z) 
parameter  (four=1.0d0,rr=1.0d0) 
integer*4  layer(mdim) 
real*  8  h(mdim) 
real*  8  rro(mlay) 

real  *8  ml  1  (mdim),ml2(mdim),mr  1  (mdim2),mr2(mdim2) 
complex*  16  lame(mlay,3) 

complex*  1 6  al  1  (mdim),al2(mdim),gl  1  (mdim),gl2(mdim) 
complex*  16  arl(mdim2),ar2(mdim2),brl(mdim2),br2(mdim2) 
complex*  16  grl(mdim2),gr2(mdim2),dll(mdim),drl(mdim) 
complex*  16  dr2(mdim),dr3(mdim),erl(mdim2),er2(mdim2) 
complex*  16  trl(mdim),tr2(mdim),qrl(mdim),qr2(mdim) 
complex*  16  prl(mdim),pr2(mdim) 
complex*  1 6  con, con  1 

all(l)-0.0d0 

gll(l)=0.0d0 

mll(l)=0.0d0 

arl(l)=0.0d0 

arl(2)=0.0d0 

brl(l)=0.0d0 

brl(2)=0.0d0 

grl(l)=0.0d0 

grl(2)-0.0d0 

mrl(l)=0.0d0 


mrl(2)=0.0d0 

erl(l)=O.OdO 

drl(l)=O.OdO 

dll(l)=O.OdO 

trl(l)=O.OdO 

qrl(l)=O.OdO 

prl(l)=O.OdO 

do  i=l,nsub 
k=2*i 
j=layer(i) 

con=^h(i)*  lame(j  ,2)/3 . 

con  1  =h(i)  *  (2 .  *  lame(j  ,2)+lame(j ,  3  ))/3 . 

all(i)=all(i)+con 

all(i+l)=con 

al2(i)=0.5*con 

ar  1  (k- 1  )=arl  (k- 1  )+con  1 

arl(k)=arl(k)+con 

arl(k+l)=conl 

arl(k+2)=con 

ar2(k-l)=0.5*conl 

ar2(k)=0.5*con 

con  1  ==lame(j  ,2)  *  h(i)/(3 .  *  rr) 

con=2.*conl 

erl(i)=erl(i)+con 

erl(i+l)=con 

er2(i)=conl 

con=conl*four 

conl=0.5*con 

trl(i)-trl(i)-i-con 

trl(i+l)=con 

tr2(i)=conl 

conl=2.*con/rr 

con=2.*conl 

qrl(i)=qrl(i)+con 

qrl(i+l)=con 

qr2(i)=conl 

con=lame(j,3)*0.5 

conl=lame(j,2)*0.5 

drl(i)=(irl(i)+con 

drl(i+l)=-con 

dr2(i)=conl 

dr3(i)=con 

dll(i)=dll(i)-conl 

dll(i+l)=conl 

con=four  *  lame(j  ,2)*0.5/rr 

prl(i)=prl(i)+con 

prl(i+l)=-con 

pr2(i)=con 

con=lame(j  ,2)/h(i) 

coni =(2 .  *  lame(j  ,2)+lame(j  ,3  ))/h(i) 

gll(i)-gll(i)+con 

gll(i+l)=con 


gl2(i)=-con 

gr  1  (k- 1  )=gr  1  (k- 1  )+con 
gr  1  (k)=gr  1  (k)+con  1 
grl(k+l)=con 
grl(k+2)=conl 
gr2(k-l)=-con 
gr2(k)=-conl 
cons=rro(j)*h(i)/3. 
mIl(i)=mll(i)+cons 
mll(i+l)=cons 
ml2(i)=0.5*cons 
mr  1  (k- 1  )=mr  1  (k- 1  )+cons 
mr  1  (k)=mr  1  (k)+cons 
mrl(k+l)=cons 
mrl(k+2)=cons 
mr2(k- 1  )=0 . 5  *  cons 
mr2(k)=0.5*cons 
con=0.5*(lame(j,2)“lame(j,3)) 
con  1=0.5*  (lame(j  ,2)+lame(j  ,3 )) 
brl(k-l)=brl(k-l)+con 
brl(k)=-conl 
brl(k+l)=-con 
br2(k-l)=conl 
br2(k)=0. 
end  do 
return 
end 

subroutine  hspace  (all,al2,gll,gl2,ml2,arl,ar2,brl,br2, 

1  gr  1  ,gr2,mr2,dr2,dr3  ,er2, 

2  tr2,qr2,pr2, 

3  oni,rock, 

4  layer, 

5  lame, vs, anu, damp, 

6  nodes,nodes2, 

7  mlay,mdim,mdim2) 
implicit  real*  8  (a-h,o-z) 

real*  8  rock 
integer*4  layer(mdim) 
real*8  vs(mlay),anu(mlay),damp(mlay) 
real*8  ml2(mdim),mr2(mdim2) 
complex*  16  lame(mlay,3) 

complex*  16  all(mdim),al2(mdim),gll(mdim),gl2(mdim) 
complex*  16  dr2(mdim),dr3(mdim) 
complex*  16  er2(mdim),tr2(mdim) 
complex*  16  qr2(mdim),pr2(mdim) 

complex*  1 6  ar  1  (mdim2),ar2(mdim2),br  1  (mdim2),br2(mdim2) 
complex*  16  grl(mdim2),gr2(mdim2) 
complex*  16  cp,cs,cpcs,coef 


calculates  half-space  contribution  (tranmitting  boundary) 


c 
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if(rock.eq.l)  then 
j  =  layer(nodes) 

cs  =  vs(j)*cdsqrt(dcmplx(1.0d0,2.0d0*damp(j))) 

cp  -  cs*dsqrt((2.0d0-2.0d0*anuG))/(1.0d0-2.0d0*anu(j))) 

cpcs  =  cp/cs 

coef  =  (0.0d0,L0d0)*lameG,2)*cs/2.0d0/om 
all  (nodes)  =  all  (nodes)  -  coef 

arl(nodes2-l)  =  arl(nodes2-l)  +  coef*'(1.0d0-2.0d0*cpcs) 
arl(nodes2)  =  arl(nodes2)  +  coef''cpcs*cpcs*(cpcs-2.0d0) 
coef  =  (0.0d0,1.0d0)*lameG,2)*om/cs 
gll  (nodes)  =  gll  (nodes)  +  coef 
grl(nodes2-l)  =  grl(nodes2-l)  +  coef 
grl(nodes2)  =  grl(nodes2)  +  coef^cpcs 
coef  =  lanieG,2)*(cpcs-2.0d0) 
brl(nodes2-l)  =  brl(nodes2-l)  -  coef 
endif 

al2(nodes)=0. 

gl2(nodes)=0. 

ml2(nodes)=0. 

ar2(nodes2-l)=0. 

ar2(nodes2)=0. 

brl(nodes2)=0. 

br2(nodes2“2)=0. 

br2(nodes2-l)=0. 

br2(nodes2)=0. 

gr2(nodes2-l)=0. 

gr2(nodes2)=0. 

mr2(nodes2“l)=0. 

mr2(nodes2)=0. 

er2(nodes)=0. 

dr2(nodes)=0. 

dr3(nodes)=0, 

tr2(nodes)=0. 

qr2(nodes)=0. 

pr2(nodes)=0. 

return 

end 

subroutine  prelove  (bll,bl2,gll,gl2,mll,ml2, 

1  oms, 

2  nodes, 

3  mdim) 

implicit  real*  8  (a-h,o-z) 
real*8  mil  (mdim), ml2(mdim) 

complex*  16  bll(mdim),bl2(100),gll(mdim),gl2(mdim) 

do  i  =  1, nodes 
bll(i)  =  oms*mll(i)  -  gll(i) 
bl2(i)  =  oms*ml2(i)  -  gl2(i) 


end  do 
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return 

end 

subroutine  love  (al,a2,bl,b2, 

1  e,v, 

2  b22,p,q,r,w, 

3  mt,n, 

4  mdim,mdim2,mdini3,*) 
c 

c  - 

c  solves  eigenvalue  problem  for  love  waves 

c  - - 

c 

implicit  real*  8  (a-h,0“z) 
integer*4  int(mdim3) 
complex*  16  c,d,ev,s,x,y,z 

complex*  1 6  a  1  (mdim),a2(mdim),b  1  (mdim)5b2(mdim) 
complex*  16  b22(mdim3),e(mdim),p(mdim2),q(mdim2),r(mdim2) 
complex*  16  w(mdim2),v(mdim,mdim) 
data  tau/1.0d~7/,eps/1.0d-9/ 

tol=cdabs(a  1  ( 1  ))*  1  .Od-20 
nl=n-l 
n2=n-2 
do  20  j=l,n 
20  b22(j)=b2(j)*2.0d0 
t=0.0d0 
ev=t 

c  search  for  root  of  determinant  by  newton  method 
do  500  nv=l,n 
nvl=nv-l 
tp=l. 
it=0 
nit=0 

ev=ev*.999 
goto  108 
105  ev=ev*(L+tau) 
tp=0. 

108  it-it+1 
if(it.lt.l00)  go  to  109 

write  (*,*) '  failure  to  converge  in  100  newton  steps' 
return  1 

109  d=0. 

if(nvLeq.0.0d0)  go  to  1 15 
do  110j=l,nvl 
c-ev-e(j) 

t=dabs(dreal(c))+dabs(dimag(c)) 
if(t.lt.tau*(dabs(dreal(e(j)))+dabs(dimag(eQ)))))  go  to  105 
110d=d-l./c 
115  x=al(l)-ev*bl(l) 

if(dabs(dreal(x))+dabs(dimag(x)).lt.tol)  x=tol 


y=-bl(l)/x 
d=d+y 
do  120  j=2,ii 

i=j-l 

z=a2(i)-ev*b2(i) 

c=z/x 

x=a  1  (j)-ev*b  1  (j)-z*c 

if(dabs(dreal(x))+dabs(dimag(x)).lt.tol)  x=tol 
y=(-b  1 0)+(b22(i)+z*y)*c)/x 
120  d=d+y 
d=-l./d 
ev=ev+d 
c=d/ev 

t=dabs(dreal(c))+dabs(dimag(c)) 

if(t.gt.eps)  go  to  108 

if(nit.eq.  1)  go  to  200 

nit=l 

tp=t 

go  to  108 
200  continue 
e(nv)=ev 

c  compute  eigenvector  by  inverse  iteration 
do  210  j=l,n 
210  w(j)=l. 

c  reduction  of  matrix 
x=al(l)-ev*bl(l) 
y=a2(l)-ev*b2(l) 
r(l)=y 

do  250  i=l,nl 
j=i+l 

c=al(j)-ev*bl(j) 

s=r(i) 

r0)=a2(j)-ev*b2(j) 
t=dabs(dreal(x))+dabs(dimag(x)) 
if(t.ge.dabs(dreal(s))+dabs(dimag(s)))  go  to  230 
p(i)=l./s 
q(i)-c 
int(i)  =  1 
z=-x*p(i) 
x=y+z*c 
y=z*r(j) 
goto  250 
230  if(t.lt.tol)  x=tol 
p(i)=l./x 

q(i)=y 

int(i)  =  0 
z=-s*p(i) 
x=c+z*y 

y=r(j) 

250  v(i,nv)=z 

if(dabs(dreal(x))+dabs(dimag(x)).lt.tol)  x=tol 
niter=0 

300  niter=niter+l 
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c  reduction  of  vector 

do  330  i=l,nl 
z=v(i,nv) 

if(int(i).eq.l)  go  to  320 
w(i+ 1  )=w(i+ 1  )+z*w(i) 
go  to  330 
320  y=w(i) 
w(i)=w(i+l) 
w(i+l)=y+z*w(i) 

330  continue 
c  back  substitution 

w(n)=w(n)/x 

w(n  1  )=(w(n  1  )-q(n  1 )  *  w(n))  *p(n  1 ) 

do  350  1=2, nl 

i=n-l 

y=w(i)-q(i)*  w(i+ 1 ) 
if(int(i).eq.l)  y=y-r(i+l)*w(i+2) 

350  w(i)=y*p(i) 
if(niter.lt.2)  go  to  300 
c  normalize  eigenvector 
cons=cdabs(w(l)) 
do  1  j=2,n 

1  if(cons.lt.cdabs(w(j)))  cons=cdabs(w(j)) 
cons=l./cons 

do  2  j=l,n 

2  w(j)=w(j)*cons 
c=w(l)*(bl(l)*w(l)+b2(l)*w(2)) 
do  420  j=2,nl 

420  c=c+w®*(b2(j-l)*w(j-l)+bl(j)*w(j)+b2(j)*w(j+l)) 
c=c+w(n)  *  (b2(n  1 )  *  w(n  1  )-i-b  1  (n)  *  w(n)) 
c=l./cdsqrt(c) 
do  440j=l,n 
440  v(j,nv)=w(j)*c 
500  continue 
do  51  i=l,n 
c=cdsqrt(e(i)) 
cr=dreal(c) 
ci=dimag(c) 

if(ci.gt.0.0d0.or.(ci.eq.0.0d0.and.cr.lt.0.0d0))  c=-c 
51  e(i)=c 

return 

end 

subroutine  prerayl  (crl,cr2,grl,gr2,mrl,mr2, 

1  oms, 

2  nodes, 

3  mdim2) 

implicit  real*  8  (a-h,o-z) 

real*  8  mrl(mdim2),mr2(mdim2) 

complex*  16  crl(mdim2),cr2(mdim2),grl(mdim2),gr2(mdim2) 
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do  i  ==  1, nodes 
kl=2*i 
k2  =  kl-l 

crl(k2)  =  grl(k2)  -  oms*nirl(k2) 
crl(kl)  =  grl(kl)  -  oms*mrl(kl) 
cr2(k2)  =  gr2(k2)  -  oms*mr2(k2) 
cr2(kl)  =  gr2(kl)  -  oms*mr2(kl) 
end  do 

return 

end 

subroutine  raylgh  (al,a3,b2,b4,cl,c3, 

1  e,v, 

2  a,v  1  ,v2,u  1  ,u2,r  1  ,r2,s  1  ,s2, 

3  mvb,nn, 

4  mdim2,mdim3,*) 
c 

c  solves  eigenvalue  problem  for  rayleigh  waves 

c  - 

c 

implicit  real*  8  (a-h,o-z) 
integer*4  mvb(mdim3) 
complex*  16  cf 

complex*  16  e(mdim2),v(mdim2,mdim2) 

complex*  16  rl(mdim2),r2(mdim2),ev,evs,ce,cd,dev,c,d 

complex*  16  a(mdim2*4),vl(mdim3),v2(mdim2),ul(mdim2),u2(mdim2) 

complex*  1 6  s  1  (mdim2),s2(mdim2) 

complex*  1 6  al  (mdim2),a3  (mdim2),b2(mdim2),b4(mdim2) 

complex*  1 6  c  1  (mdim2),c3  (mdim2) 

data  epsl/l.e-5/,eps2/l.e-10/ 


nl=nn-l 
n2=nn-2 
do  10j=l,n2,2 
mvb(j)^j+3*nn 
H+i 

10  mvb(i)=i+2*nn 
mvb(nl)=nl+nn 
xn=0. 

do  15  j=l,nn 
x=cl(j)/al(j) 

15  xn=xn+dabs(x) 
xn=dsqrt(xn/nn) 
c=dcmplx(xn,xn+xn) 
do  20  j=l,nn 

sl(j)=l. 

s2(j)=c 

vl(j)=sl0) 

v20)=s2a) 

rl(j)=0. 

20  r2G)=0. 


ev=xn/nn 

isw=0 

kc=0 

mc=0 

400  kc=kc+l 
do  50n=l,nn 
ul(n)=cl(n)*vl(n) 

50  u2(n)=al{n)*v2(n) 
do  55  n=3,nn 
l=n-2 

u  1  (l)=u  1  (l)+c3  (1)  *  V 1  (n) 
ul  (n)=u  1  (n)+c3(l)*  V 1  (1) 
u2(l)=u2(l)+a3(l)*v2(n) 

55  u2(n)-u2(n)+a3(l)*v2(l) 
if(isw.eq.l)  go  to  85 
isw=0 
cd=l. 

x=cdabs(ev)/2. 
ev=dcmplx(x,x+x) 
go  to  87 

85  ev=conjg(ev) 
isw=2 
cd=l. 

87  ik=0 

do200it=l,100 

evs=ev*ev 

do  90  n=l,nn 

a(n)=evs  *  a  1  (n)+c  1  (n) 

a(n+nn)=ev*b2(n) 

a(n+2  *  nn)=evs  *  a3  (n)+c3  (n) 

a(n+3  *nn)=ev*b4(n) 

90  v2(n)=(ul(n)-u2(n)*ev)*cd 
do  150  n=l,nl 
i=n 

j=n+nn 
m=mvb(n) 
do  120  l=j,m,im 
c=a(l)/a(n) 
i=i+l 

ji=i 

do  1 10  k=l,m,nn 
a(ji)=a(ji)-c*a(k) 

110  ji=ji+im 
a(l)=c 

120  v2(i)-v2(i)-c*v2(n) 

150  v2(n)=v2(n)/a(n) 
n=nn 

if(cdabs(a(n)).eq.0.0d0)  then 

dev=dev*0.5 

ev=ev-dev 

goto  200 

endif 

v2(n)=v2(n)/a(n) 
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do  160  k=l,nl 

n=n“l 

i=n 

j=n+nn 
m=invb(n) 
do  160  l=j,m,nn 
i=i+l 

160  v2(n)=v2(n)-a(l)*v2(i) 
ce=0. 

do  170  n=l,nn 
V 1  (n)=(v2(n)- V 1  (n)*  cd)/ev 
ce=ce-u  1  (n)*  V 1  (n)+u2(n)*  v2(n) 
ul(n)=cl(ii)*vl(n) 

170  u2(n)=al(n)*v2(n) 
ce=ce*cd 
do  180  n=3,nn 
l=n-2 

ul(l)=ul(l)+c3(l)*vl(n) 
ul(n)=ul(n)+c3(l)*vl(l) 
u2(l)=u2(l)+a3  (1)  *  v2(n) 

180  u2(n)=u2(n)+a3(l)*v2(l) 
cd=0. 

do  190  n=l,nn 

190  cd=cd-ul(n)*vl(n)+u2(n)*v2(n) 
dev=ce/cd 

if(it.gt.l5)  dev=dev*0.5 

ev=ev+dev 

cf=cd 

cd=cdsqrt(2./cd) 

c=dev/ev 

x=dabs(dreal(c))+dabs(dimag(c)) 
if(ik.eq.l.and.x.lteps2)  go  to  300 
if(x.lt.epsl)  ik=l 
200  continue 

write  (*,*) '  failure  to  converge  in  100  iteration  steps' 
return  1 
300  e(kc)=ev 
icr=0 

x=dabs(dreal(ev)) 

y=dabs(diniag(ev)) 

z==x+y 

if(y/z.lt.eps2)  icr=l 
if(x/z.lt.eps2)  icr=2 
if(icr.eq.l)  e(kc)=dreal(ev) 
if(icr.eq.2)  e(kc)=diniag(ev)*(0.,  1 .) 
do  310  n=l,nn 
310  v(n,kc)=v2(n)*cd 

c  if  the  eigenvalue  is  real  choose  the  sign  such  that 
c  the  group  velocity  becomes  positive 
if(icr.ne.l)goto316 
x=(vl(l)*vl(l)+vl(2)*vl(2))*ev 
if(x.gt.0.0d0)  go  to  3 16 
e(kc)=-e(kc) 


do  3 15  n=2,nn,2 
315  v(n,kc)=-v(n,kc) 

c  check  the  special  orthogonality  of  the  newly  found  eigenvector 
c  and  the  sum  of  the  previously  found  eigenvector 
316ce=0. 
do  330  n=l,nn 

330  ce=ce-rl(n)*ul(n)+r2(n)*u2(n) 
ce=ce*cd 

x=dabs(dreal(ce))+dabs(dimag(ce)) 
if(x.lt.epsl)  go  to  335 

write  (*,*)  ’  failure  to  find  an  orthogonal  eigenvector’ 
return  1 
335  c=l./e(kc) 
do  340  n=l,nn 
r  1  (n)=r  1  (n)+v(n,kc)*  c 
340  r2(n)=r2(n)+v(n,kc) 

c  find  a  starting  vector  orthogonal  to  all  eigenvectors  found 
c  for  the  iteration  toward  the  next  eigenvector 
c=0. 
d=0. 

do  342  n=2,nn,2 
m=n-l 

c=c-u  1  (m)*s  1  (m)+u2(n)*  s2(n) 

342  d=d-ul(n)*sl(n)+u2(m)*s2(m) 
ce=2./cf 

c=c*ce 

d=d*ce 

do  343  n=2,nn,2 
m=n-l 

s  1  (m)=s  1  (m)-c*  V 1  (m) 

s  1  (n)=s  1  (n)-d*  V 1  (n) 

s2(m)=s2(m)-d*v2(m) 

s2(n)=s2(n)-c*v2(n) 

vl(m)=sl(m) 

vl(n)=sl(n) 

v2(m)=s2(m) 

343  v2(n)=s2(n) 
if(icr+isw.eq.O)  isw=l 
mc=mc+l 

if(mc.ge.nn)  go  to  350 
go  to  400 

c  form  vectors  for  displacement  expansion 
350  do610n=l,nn 
c=-l. 

x=dimag(e(n)) 
if(x.le.0.)  go  to  600 
e(n)=-e(n) 
c=-c 

600  do  610  j=2,nn,2 
610  vQ,n)=vG,n)*c 


return 

end 


SUBROUTINE  HANKEL(ZZ,H0,H1,IND) 


C - 

C  Computation  of  zero  and  first  order  Hankel  functions 
C 

C  ZZ  -  complex  argument,  -3.1415...LE.ARG(Z)  .LE.  3.1415... 

C  HO  -  Hankel  ftinction  of  IND'th  kind  (IND=1,2)  and  zero  order 
CHI  -  Hankel  function  of  IND'th  kind  (IND=1,2)  and  first  order 
C 

C  Written  by  G.  Waas  and  E.  Kausel 


IMPLICIT  REAL*8  (A-H,0-Z) 

COMPLEX*  16  Z,H0,H1,C,A,E,E2,ZH,P 

COMPLEX*  16  G0,G1,ZZ 

Z=ZZ 

SG=LD0 

IF(IND.EQ.1)SG=-LD0 

X=DREAL(Z) 

Y=DIMAG(Z)*SG 

R=DSQRT(X*X+Y*Y) 

PHI=DATAN2(Y,X) 

IF(R.LE.10.D0)GO  TO  6 
C  Use  asymptotic  expansion  for  R  >  10 
J=2.D0*R 
11=3 

H0=0.D0 

H1=0.D0 

IF(PHI.LT.1.57DO)GO  TO  2 

11=0 

X=-X 

Z=-Z 

SG=-SG 

PHI=3 . 1 4 1 592653 589793D0-PHI 

1  Y=-Y 
PHI=-PHI 
SG=-SG 
11=11+1 

2  G0=2.D0*H0 
G1=2.D0*H1 

C=DCMPLX(0.D0,SG*0. 1 25D0)/Z 

K=2 

P=C*C 

A=4.5D0*P 

P=7.5D0*P 

H0=LD0+C+A 

H1=1.D0-3.D0*C-P 

3  I=4*K 
K=K+1 
DI=I 
DK=K 

A=A*C*(DI+LD0/DK) 

P=P*C*(DI-3.D0/DK) 

H0=H0+A 


H1=H1-P 

IF(DABS(DREAL(P))+DABS(DIMAG(P)).GT.1.D-16.AND.K.LT.J)G0  TO  3 
C=H0/H  1  *DCMPLX(0.D0,-SG) 

AR=(O.785398163397448D0-X-PHI/2.D0)*SG 

E=0.797884560802865D0/DSQRT(R)*DEXP(Y)*DCMPLX(DCOS(AR),DSIN(AR)) 

IF(X.NE.O.OR.PHI.GT.O.)GO  TO  4 

E=DCMPLX(O.DO,DIMAG(E)) 

4  H0=H0*E 

H1=H1*E*DCMPLX(0.D0,SG) 

HO=HO+GO 
H1=G1+H1 
GO  TO  (1,5, 9), II 

5  H1=-H1 
GO  TO  8 

C  Use  ascending  series  for  R  <=  10 

6  ZH=Z/2.D0 
C=-ZH*ZH 

E=DCMPLX(0.D0,.3 1 8309886183791D0)*SG 
E2=E*2.D0 

A=l.D0-E2*(0.577215664901533D0+DLOG(R/2.D0))+PHl*0.636619772367582 

$D0 

P=1.D0 

K=1 

H0=A 

H1=A+E*(1.D0-1.D0/C) 

7  A=A+E2/K 
P=P*C 
H0=H0+A*P 
K=K+1 
P=P/(K*K) 

H1=H1+(A*K+E)*P 

IF(DABS(DREAL(P))+DABS(DIMAG(P)).GT.1.D-16)  GO  TO  7 
H1=H1*ZH 

IF(X.NE.O.OR.PHI.GT.O.)GO  TO  8 
H0=DCMPLX(0.D0,DIMAG(H0)) 

H1=DREAL(H1) 

8  C=H0/H1 

9  RETURN 
END 

SUBROUTINE  BESSEL  (ZZ,J0,JI) 

C - 

C  Computation  of  zero  and  first  order  Bessel  functions 
C 

C  ZZ  -  complex  argument,  -3.1415...LE.ARG(Z)  .LE.  3.1415... 

C 

C  Written  by  E.  Kausel 

C - 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMPLEX*  16  Z,ZZ,C,P,A,AO,A1,BO,B1,JO,J1 
Z  =  ZZ 

X  =  DREAL(Z) 

Y  =  DIMAG(Z) 


o  o  o 


R  =  DSQRT(X*X  +  Y*Y) 
IF  (R.LE.10.D0)  GO  TO  30 


Use  asymptotic  expansion  for  R  >  10 

PHI  =  DATAN2  (Y,X) 

L  =  0 

IF  (X.GE.0.D0)  GO  TO  10 
L=  1 
Z  =  -Z 

PHI  =  PHI  -  3.14159265358793D0  *  DSIGN  (l.D0,PHI) 

X  =  -X 
Y  =  -Y 
10  J  =  2.D0*R 

C  =  (O.DO,-0.125DO)/Z 
K  =  2 
P  =  C*C 
A  =  4.5D0*P 
P  =  7.5D0*P 
A0  =  C 

B0  =  1.D0  +  A 
A1  =-3.D0*C 
B1  =  l.DO-P 
20  I  =  4*K 
K  =  K+1 
DI  =  I 
DK  =  K 

A  =  A*C  *  (DI  +  1.D0/DK) 

P  =  P*C  *  (DI  -  3.D0/DK) 

AO  =  AO  +  A 
A1  =A1  -P 
I  =  4*K 
K  =  K+1 
DI  =  I 
DK  =  K 

A  =  A*C*  (DI+  l.DO/DK) 

P  =  P*C  *  (DI  -  3.D0/DK) 

BO  =  BO  +  A 
B1  =B1  -P 

IF  (DABS(DREAL(P))+DABS(DIMAG(P)).GT.1.D-16  AND.  K.LT.J)  GO  TO  20 
CC  =  DCOS  (X) 

S  =  DSIN  (X) 

SS  =  S  +  CC 
CC  =  S  -  CC 
TH  =  DTANH  (Y) 

CH  =  0.56418958354776D0  *  DCOSH  (Y)/DSQRT  (R) 

S  =  PHI/2.D0 

A  =  CH  *  DCMPLX  (DCOS(S),-DSIN(S)) 

JO  =  BO  *  DCMPLX  (SS,-CC*TH) 

JO  =  JO  +  A0*DCMPLX  (-SS*TH,CC) 

J1  =  B1  *  DCMPLX  (CC,SS*TH) 

J1  =  J1  -  A1*DCMPLX  (CC*TH,SS) 

IF  (L.EQ.1)  J1=-J1 


o  o  o 
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C  =  JO/Jl 
JO  =  JO*A 
J1  =J1*A 
RETURN 


Use  ascending  series  for  R  <=  10 

30  A  =  Z/2.D0 
C  =  -A*A 
J0  =  l.DO 
J1  =  l.DO 
P=  l.DO 
K=  1 

40  P  =  P*C 
JO  =  JO+P 
K  =  K+1 
DK  =  K 
P  =  P/DK 
J1  =  J1  +P 
P  =  P/DK 

IF  (  DABS  (DREAL(P))  +  DABS  (DIMAG(P))  .GT.  l.D-16)  GO  TO  40 
J1  =  J1*A 
c  C  =  JO/Jl 
RETURN 
END 


